{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "f9d691d4",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import integrate\n",
    "import math\n",
    "import numpy as np, matplotlib.pyplot as plt\n",
    "from sklearn.metrics import mean_absolute_error as MAE\n",
    "import time\n",
    "from tqdm import trange\n",
    "from environment import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f85e98ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "#import data\n",
    "import pandas as pd\n",
    "tcga_all=pd.read_csv('tcga_all.csv')\n",
    "\n",
    "#init environment\n",
    "env= environment()\n",
    "\n",
    "# initialize r_l\n",
    "num_cases=len(tcga_all)\n",
    "r_l=np.ones(num_cases)*0.1\n",
    "\n",
    "# target tcd8+ cells proportion\n",
    "tcd8_target=tcga_all['tcd8_pop'][:num_cases]\n",
    "\n",
    "#init running time\n",
    "tcga_all['time']=int(600)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "872b08ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def opt_fix_para(r_l,jc_scale,h,q,num_cases,tcd8_target):\n",
    "    cd8_l=[]\n",
    "    time_l=[]\n",
    "    for i in trange(num_cases):\n",
    "        mu=tcga_all.loc[i,'mu']\n",
    "        r=r_l[i]\n",
    "        jc=tcga_all.loc[i,'mhc_norm']\n",
    "        jc_=(1-np.exp(-jc*jc_scale))/(1+np.exp(-jc*jc_scale))\n",
    "        R_L_prop=tcga_all.loc[i,'Treg_pro']\n",
    "        tumor=tcga_all.loc[i,'tumor_size']\n",
    "        a=tcga_all.loc[i,'a']\n",
    "        ar_ratio=tcga_all.loc[i,'ar_ratio']\n",
    "        t1=int(tcga_all.loc[i,'time'])\n",
    "        \n",
    "        x_=env.fit_pre_treatment(mu,r,jc_,R_L_prop,h,q,a,ar_ratio,t1)\n",
    "        j=cal_time(x_,tumor)\n",
    "        cd8_pred=min(x_[j][2]/(x_[j][0]+x_[j][1]),0.1)\n",
    "        cd8_l.append(cd8_pred)\n",
    "#         print(int(j/24)+10,x_[j])\n",
    "        time_l.append(int(j/24)+10)\n",
    "    return MAE(cd8_l,tcd8_target),cd8_l,time_l\n",
    "def cal_time(x_,tumor):\n",
    "    for i in range(len(x_.T[0])):\n",
    "        if x_.T[0][i]+x_.T[1][i]>tumor:\n",
    "            break\n",
    "    return i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac8e872b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#maximization step: given hidden parameter r, maximize the j1_scale, h, q\n",
    "test_range=[]\n",
    "for j1_scale in [2,4,6,8]:\n",
    "    for h in range(1,4):\n",
    "        for q in range(1,4):\n",
    "            test_range.append((j1_scale,h*1e-2,q*1e-2))\n",
    "\n",
    "\n",
    "# first round, the minimum mae 0.009506 with (j1_scale,h,q)= 6 0.01 0.03\n",
    "for j1_scale,h,q in test_range:\n",
    "    mae, cd8_l,time_l = opt_fix_para(r_l,j1_scale,h,q,num_cases,tcd8_target)\n",
    "    tcga_all['time']=time_l\n",
    "    print(j1_scale,h,q,mae)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "3f4f85c4",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAACJ6klEQVR4nOzdd3gUVRfH8e/uZtMbqZCQQOgQmgTpoCKCWAAFRLErKGKh2EB8rSCKBQuCgIgVRZqAooAiUqX3DgHSSSO9bJv3j4FApKVsMpvkfJ6Hh93Z2ZkTSvaXO3PP1SmKoiCEEEIIoRG91gUIIYQQomaTMCKEEEIITUkYEUIIIYSmJIwIIYQQQlMSRoQQQgihKQkjQgghhNCUhBEhhBBCaErCiBBCCCE05aR1ASVhs9lISEjAy8sLnU6ndTlCCCGEKAFFUcjOziYkJAS9/srjH1UijCQkJBAWFqZ1GUIIIYQog9jYWOrWrXvF16tEGPHy8gLUL8bb21vjaoQQQghREllZWYSFhRV9jl9JlQgj5y/NeHt7SxgRQgghqphr3WIhN7AKIYQQQlMSRoQQQgihKQkjQgghhNBUlbhnpCQURcFisWC1WrUupUIYDAacnJxkarMQQohqp1qEEZPJRGJiInl5eVqXUqHc3d2pU6cOzs7OWpcihBBC2E2VDyM2m42TJ09iMBgICQnB2dm52o0eKIqCyWQiJSWFkydP0rhx46s2jxFCCCGqkiofRkwmEzabjbCwMNzd3bUup8K4ublhNBo5ffo0JpMJV1dXrUsSQggh7KLa/HhdE0YKasLXKIQQouaRTzchhBBCaKrUYWTdunXceeedhISEoNPp+OWXX675nn/++YeoqChcXV1p0KABX3zxRVlqFUIIIUQ1VOowkpubS5s2bZg2bVqJ9j958iS33XYb3bt3Z9euXbzyyis899xzLFq0qNTFCiGEEKL6KfUNrH379qVv374l3v+LL74gPDycjz/+GIDmzZuzfft2PvjgAwYOHFja01cbkydPZvHixRw+fBg3Nze6dOnCe++9R9OmTbUuTQghhKhUFT6bZvPmzfTu3bvYtj59+jBnzhzMZjNGo/GS9xQWFlJYWFj0PCsrq6LLrHT//PMPTz/9NNdffz0Wi4UJEybQu3dvDh48iIeHh9blCSGEKAPFZgOr9fK/WyyXbrfawGZFsVr/8/yi362W0m23WVEs1stvt9pQrBa4zHafuwbgFhmpyZ9bhYeRpKQkgoODi20LDg7GYrGQmppKnTp1LnnP5MmTefPNN8t0PkVRyDdr04XVzWgocY+TP/74o9jzuXPnEhQUxI4dO+jRo0dFlCeEENWCYrWimM0oJlOxXzaTCcV00Xbz5V7/zz4X7XvF1/9zLpv5on3M5mLhAkXR+o+nzNyua1t9wwhcunSwcu4v60of3OPHj2fs2LFFz7OysggLCyvRufLNVlq8trKMlZbPwbf64O5ctj/SzMxMAPz8/OxZkhBCVCrFYsGanY0tO/vc7zlYs7OwZedgy8nGmqW+ZsvLu3YIuNyHv8mkfuhXVQYDOoNB/V2vv/R3J6fLb7fH64bzz4ufH70OMzb0EfU0+2Op8DBSu3ZtkpKSim1LTk7GyckJf3//y77HxcUFFxeXii7NYSiKwtixY+nWrRstW7bUuhwhRA2lKAq23Dxs2VlXDhTZWVizcy567aLfc3JQNFiWQ+fsfJlfRnTOzuiNl3vt3OvGc/tcvP2S/S+zT7F9jZf50L9K2LADs81MnjmPPHMeueZcci25Rc/zLOe2mXPJs1xm27nn59+bZ8kj35IPwFSPdvRCm8+gCg8jnTt3Zvny5cW2rVq1ivbt21/2fpHycjMaOPhWH7sft6TnLotnnnmGvXv3smHDBjtXJISoSWyFhWooyMpWRyEuDgtZ2VhzsosHiqwsrDkXBYucHLDZ7FKLzs0Ng5cXei+vi373RO/ljcHLE527+1VCwFUCwEX7652NYDQ69BIgiqJQaC0ktzCrROHgSiEi35Jf9NhsM1dIrbnm3Ao5bkmUOozk5ORw/PjxoucnT55k9+7d+Pn5ER4ezvjx44mPj+fbb78FYMSIEUybNo2xY8cyfPhwNm/ezJw5c/jxxx/t91VcRKfTlflSiRaeffZZli1bxrp166hbt67W5QghHIhisWBOOoM5LhZzXBym2DgsaalXDBSKyWSfEzs5XSZIXBoo9F7e6L08L93H0xNdBfywqRWbYiOrMIu0gjTSC9JJy08jrSCNtPw0MgszybWcCwzm/MuOVFiVirms5Kx3xt3ojofRAzcnNzyMHngYPXB3csfd6I67k/uFbeeen9//cvu5GLS7IlHqT+3t27dz0003FT0/f2/Hww8/zNdff01iYiIxMTFFr0dERLBixQrGjBnD559/TkhICJ9++mmNntYLalp+9tlnWbJkCWvXriUiIkLrkoQQGrBmZmKKjcMcF4spNhZzbJwaPOLiMCckgMVS6mPqPT2vECQuHyjU3y/so3N1dejRBnuw2CycLTirBoz89KJwcbnAcbbgLBal9H8P/+Xm5Fb0wX9xCDgfEC4OC0WvXW7bufcaDdUn8JU6jNx4441FN6Beztdff33JthtuuIGdO3eW9lTV2tNPP828efNYunQpXl5eRffV+Pj44ObmpnF1Qgh7UUwmzAkJFwJHXBzm2DhMcWrwsGVnX/X9OqMRY2goxrp1MYbVxRgUdNkgcT5o6D081HsWaqBCa6EaIvLPBYqLA8ZFgSO9IJ2MwgwUSjfzxdvZGz9XP/zd/PF39cffzR9fF98LgcLpomBxcdhwcsfNyQ2Dvmb+vZRE1bmeUc3MmDEDUMPdxebOncsjjzxS+QUJIcpEURSsaWnqqEZc/LkRjjjMsbGY4uOwJCZdc7qnITAA59C6GMPCcA6ri7FuGMa6oTiHheEUFFRjw4WiKOSac4uFiItHLP4bOEp7z4Nep8fXxbdYuPBz9St67O/qj5/bueeu/tVqJMLRSBjRyNVGl4QQjsWWn485Pr7oMoop7lzwiFVHOpT8/Ku+X+fqeiFkhNXFuW4Yxrp11W2hoejd3SvpK3EMeeY8EnISLntJ5OLn6QXpFFoLr33Aixj1xktGL4oFjIue+7r4ymiFg5AwIoSo8RSbDUty8rlwcT5kXAge1pTUqx9Ap8Opdm2c66qjG+dHNdTAEYbB37/a34PxX4qikFaQxsnMk0RnRBOdGa0+zozmTN6ZUh3Lzcnt8qMV55+fDx9u/ngZvWrcn3V1IGFECFEjWHNyzs1IOX+T6LnLKXHqL8V89emSek9PjOFhl1xOcQ6ri1NICHpn50r6ShyL1WYlPie+KGhcHDqyTVe+H8bL2YsAt4BrXh7xc/XD3VizRo5qIgkjQohqyZqdTe6mzeSsX0fuho1Y/tN88RJOThhDQnCuG3rhckpYGMZQ9XKK3senRv/EXWAp4FTWqQuhIyOak1knOZ15GpPt8lOKdegI9QylgW8DGvg0IMInouh3HxefSv4KhCOTMCKEqBYURaHw8GFy1q0nd/168nbvvmRarKFWLXVU45LLKWEYawejc5JviRkFGcVGN84/TshJuOLsE2e9M/V96tPA50LoiPCJoL5PfU17V4iqQ/7nCSGqLGtmJrmbNxcFEEtKSrHXnSMi8OzRHY/uPXBr2waDp6dGlToWm2IjKTepeOjIiOZU1inSC9Kv+D4fF59LRjgifCII8QiRG0FFuUgYEUJUGYrNRsGhQ+SuX0/OuvXk79lTbNE0nZsbHp064dG9G549euBcw7sam61mTmedLhY6Tmae5FTWqaL1SC6njkedYmHj/GM/V78afalKVBwJI0IIh2bNyCB30yZy1q0nZ8MGrKnFZ7Y4N2yIZ/fuePbojlv79jXyRtJsU3axsHH+97jsuCu2InfSO1HPqx4NfBtQ37t+0X0d9b3ryw2jotJJGBFCOBTFZqPgwEH1xtN168nfu7fY4m16d3fcO3dWA0j3bhhDQzWstnJlFGRw+Oxh9ebRzJNFwSMlP+WK7/EwelwyytHApwGhXqEY9dLESwAFWZByBAKbgKs2NxZLGBFCaM5y9iy5GzaSu2E9ORs2Yk1LK/a6S+NGeHTvgWeP7ri3a4euhox+WGwW9qXuY0P8BjbFb+JA2oEr3kQa5BZUFDgifCKKRjoC3QLl0opQmXLV0JF8CFIOQfJhSDkMmbHq6/fNh6a3alKahBEHMXnyZF555RVGjRrFxx9/rHU5QlQoxWql4MAB9dLL+nUU7N1XrGW63sMDjy6d8ejeHc9u3TCGhGhYbeVKyk1iY/xGNiZs5N+Ef8k2F+/VEe4Vftmpsl7OXhpVLByOOR9Sj54LG4fU8JF8CDJi4Err8XjWBlNOpZZ5MQkjDmDbtm3MmjWL1q1ba12KEBXGkp5O7oYN6syXjRuxnj1b7HWXJk2KZr64X9e2xox+FFoL2ZG0g40JG9kYv5ETmSeKve7j4kPnOp3pGtqVLiFdCHIP0qhS4XAshZB6TB3dSD507veDcPYUKLbLv8cjEAKbQVDz4r+7+1Vq6f8lYURjOTk53H///cyePZuJEydqXY4QdqNYrRTs23du9GM9Bfv3Fx/98PTEo0uXcwGkO8bgYA2rrTyKonAq6xQb4zeyIWEDO5J2UGAtKHpdr9PTKqAVXUO70jWkK5H+kTJttqazmiHthBo0Lg4eaSfgCjco41YLApurYePi4OERULm1l1D1CyOKAuY8bc5tdIdSXpt9+umnuf322+nVq5eEEVHlWVJTydmwgdzzox+ZmcVed2nW7MLMl7Zt0Rlrxg2UOaYctiRtUS+/xG8kITeh2OtBbkFq+AjtSqc6naQ7aU1ltcDZkxcuq5y/ryPtONiusFyBiw8ENTsXNlqce9wcPINK/XmkpeoXRsx58I5G15dfSQBnjxLv/tNPP7Fz5062bdtWgUUJUXEUi4X8vXvJWb+e3HXrKThwoNjrei8vPLp2xbN7dzy6dcMYXDMuMdgUG4fTDxfd+7EneQ8W5UI3WKPeSFRwFF1D1ADSyLeR3GRak9is6qWU86Mc50c6Uo+C9fKt9XH2PBc4mhUf8fCqU6VCx5VUvzBSRcTGxjJq1ChWrVqFq6ur1uUIUWKWlBRy1m9Qp95u2oztv6MfLZrjeW7mi1ubNjWmxXp6QTqbEjaxMX4jmxI2XdLJtJ53vaLw0T64vfTyqAlsNsiMuehG0nO/pxyFKzWdM7pDYNNzgaPZhd99wqpF6LiS6vddwuiujlBode4S2rFjB8nJyURFRRVts1qtrFu3jmnTplFYWIjBINeJhfYUi4X83buL7v0oPHSo2Ot6Hx88u3ZRp95264pTYKBGlVYui83C3pS9bIjfwMaEjRxKO1Rs2q27kzsd6nSgW0g3uoR2IcwrTMNqRYVSFMiKLz7KkXxInUZrzr38ewwual+P/97X4VsP9PrKrd8BVL8wotOV6lKJVm6++Wb27dtXbNujjz5Ks2bNePnllyWICM2ZTp8mZdrn5Kxdiy27+PRS18hIPHp0x7N7D9xat6oxox+JOYlFs17+TfyXHHPxqZDN/JrRJaQL3UK70TawLUZDzbgnpsbJSYHTG+H0JkjYqYaOwqzL72twBv/Gl15eqVUf5MbkIjXjO4gD8vLyomXLlsW2eXh44O/vf8l2ISqTLS+P1C9mkj53LopZvWnO4OODR7du6syXbt1w8vfXuMrKUWApYMeZHWrTsYRNRGdGF3vd18WXziGd6RqiTrsNdK8Zo0I1Tma8GjxOb1R/pR69dB+9E/g3+s+02Rbg1wAM8lF7LfInJIQA1Cmn2b//zpkp72NJSgLAo2tXAp5+Grc2rdHVgNE6RVE4mXWyaNbL9jPbKbQWFr2u1+lpHdCarqFd6RbajeZ+zWXabXWjKOrNpReHj7OnLt0vKBLqd4Wwjmro8G8ETjWjN05FkDDiQNauXat1CaKGKjhylDOTJpG3dSsAxtBQgsePw/Pmm6v9LI9sUzZbErcUXX5JzE0s9nqwe3BRz4+OdTrKtNvqRlHUxmHng8fpTer9HxfT6aF2a6jfDep1gfDOmjcJq24kjAhRg1mzskj5bBpn580DqxWdiwv+TwzH//HH0VfTWV42xcah9ENFox97UvYUW9nWWe+sTrs9F0Aa+jas9oGsRrHZ1OZhF4eP3P8sNKg3Qmg7NXjU6wZhHcDVW5t6awgJI0LUQIrNRuaSJSR/+BHWdHUKqlfv3gS//FK1XAU3vSC9qOfH5oTNl0y7re9dvyh8tK/dHjcnN40qFXZntUDS3gvB4/QmKMgovo+TK9S9/lz46Ko+dpap15VJwogQNUz+3r0kTZxEwd69ADg3aEDtVyfg0aWLxpXZX2x2LF/t/4qlx5divqiDpbuTO53qdCpa76WuV10NqxR2ZTGpM1xOb4RTGyF2y6ULwBk9ILyjGjzqdVVHQZxctKlXABJGhKgxLGlpJH/0EZmLFgPqyrgBTz+N3wP3V7tF6U5knODLfV/y+8nfiy7BNK3VlO51u9MlpItMu61OzPkQt00d8Ti1AeK2X9pQzMUH6nW+ED7qtAb5+3coEkaEqOYUi4Wz834k5bPPivqF+PTvT+DzYzEGVa/27AfSDvDl3i/5M+bPom1dQ7oyvPVwooKjrvJOUWUUZqujHac3qSMf8TsuXbfF3f/C/R71ukBwpPT0cHASRoSoxnK3bOXMxIkUHjsGqK3aa7/6P9zbXadxZfa148wOZu+bzcb4jUXbbg6/meGthhMZEKlhZaLc8s9CzL/qqMfpTZC459KVar3qnBv16KLOeAloUq1bp1dHEkaEqIbMSUkkT5lC1orfAbVpWeCYMfgOHlRt+oUoisKmhE3M2juLnck7ATDoDPSN6MvjLR+nUa1GGlcoyiQnBWLOjXqc3gRn9sNFbfYBtWV6va5qn496XaBWhISPKk7CiBDViM1kIn3u16R+8QVKfj7o9fgOuYfA557DqVYtrcuzC5ti4++Yv5m9bzYH0tRVgo16I/0b9eexlo/JGjBVTVbCueBxle6m/o0vjHrU6wI+csNxdSNhRIhqInvtWs5Mnoz5dAwAbu3aUfvVCbi2aKFxZfZhsVn449QfzNk3h+MZxwFwNbgyqMkgHol8hGCPYI0rFCViyoWjf8DxNXB6w9W7m9brAuFdwEv+bqs7CSNCVHGm06c5M/ldcs518HUKDCTopRfxvuOOatGsy2Q1sezEMr7a/xWx2bEAeBo9ua/ZfTzQ4gH8XKUTpsOzmODEX7BvIRxZAea8C69Jd1OBhBFNxcfH8/LLL/P777+Tn59PkyZNmDNnDlFRcte/uDZbXh6pM2eR/tVX6oJ2Tk74PfwQAU+NxODp+CtXX0u+JZ9FRxcx98BckvOSAXVhugdbPMi9ze7F21k6Yjo0mxVOrVcDyKFlUJB54bVa9aF5P4jooa7tIt1NazwJIxo5e/YsXbt25aabbuL3338nKCiIEydO4Ovrq3VpwsEpikL2H39w5r0pxRa0C57wCi4NGmhcXfllm7KZf2Q+3x38rqhTapBbEA9HPsygJoNwN0pnTIelKGqfj/0L4cASyDlz4TXP2tDybmg5SG0yVg1G7YT9SBjRyHvvvUdYWBhz584t2la/fn3tChJVQsHRo5yZ9A55W7YA1WtBu7MFZ/n+0Pf8eOhHss1qP5RQz1Aea/kYAxoNwNlQvRqzVStnDqgjIPsXQcbpC9tdfaFFf2g1SJ39Ir0+xBVUuzCiKAr5/+2+V0ncnNxK/IGwbNky+vTpw+DBg/nnn38IDQ1l5MiRDB8+vIKrFFWRNSuLlGnTOPvDRQvaDR+O/7Cqv6Bdcl4y3xz4hgVHFxT9323g04BhrYbRN6IvTvpq922qekiPhn2L1ACScujCdqMHNLsdWg6Ehj3BSUKkuLZq978835JPx3kdNTn3lqFbSjyEHB0dzYwZMxg7diyvvPIKW7du5bnnnsPFxYWHHnqogisVVYW6oN0vJH/0Eda0NAC8bulF0MvjcK5btRe0i8uOY+7+uSw5vqRo3Zjmfs0Z3no4N4ffjF6n17hCcYmsRDiwWB0FSdh5YbvBGRr3VgNIk1tlkTlRatUujFQVNpuN9u3b88477wBw3XXXceDAAWbMmCFhRACQv28fSRMnUrDnwoJ2wRNewbNrV40rK5/ojGjm7J/Db9G/Fa0b0y6oHcNbD6drSNcqf7mp2slLh4NL1RGQUxsoakCm00PEDeolmGZ3gJuvllWKKq7ahRE3Jze2DN2i2blLqk6dOrT4T/+H5s2bs2jRInuXJaoYS3o6KVOnkrFwESgKend3dUG7Bx+o0gvaHUo7xOx9s/nz9J8o5z7QuoR0YXir4bSv3V7j6kQxhdlweIV6I+qJNWCzXHgtrJMaQFr0B8/qtbaR0E61CyM6na5K3G3ftWtXjhw5Umzb0aNHqVevnkYVCa0pFgtnf5pPyqefYsvKAsCnfz8Cn3++Si9otyt5F7P2zmJD/IaibT3DejK89XBaBrTUsDJRjLkAjq9WL8EcXVl85dvardRZMC3vBt9w7WoU1Va1CyNVxZgxY+jSpQvvvPMO99xzD1u3bmXWrFnMmjVL69KEBvK2bSPp7YkUHlVbYbs0b07t/72Ke7t2GldWNoqisDlxM7P3zmb7me0A6HV6bq1/K8NaDaNxrcYaVygAsFrg5D/qJZhDy6Ew68Jrfg3VEZCWgyCwiXY1ihpBwohGrr/+epYsWcL48eN56623iIiI4OOPP+b+++/XujRRicxnzpA85X2yfvsNOL+g3Wh8Bw+ukgva2RQba2PXMnvvbPan7QfASe9E/4bqujHh3vJTteZsNojbqo6AHPwFclMuvOYdCpF3qSGkTlvpBSIqjYQRDd1xxx3ccccdWpchNGAzmUj/+ht1Qbu8PNDp1AXtRo2qkgvaWW1WVp5ayex9s4utGzOwyUAeiXyE2h61Na6whlMUSNqrBpADSyAz9sJr7v7QYoAaQMI6gV5mMYnKJ2FEiEqW888/nHlnMqbTanMot+uuo/b/Xq2SC9qZrWaWRy9nzr45xGSrC/R5GD3UdWOaP4C/m7/GFdZwqcfVm1D3LYS0Yxe2O3tB8zvUSzANbgCDUbsahUDCiBCVxhQToy5o9/ffABgCAwh+4QW8+/WrctNZ8y35LD62mLn753ImT2357eviywPNH+C+5vfJujFayoyD/YvVEJK458J2gws06aOOgDTuDcaSz/4ToqJJGBGigtny80mdNYv0OV+hmEzqgnYPPUTAyKcweHpqXV6p5Jhy+OnIT8XWjQl0C+ThyIcZ3GRwlZjJVi3lpqqXX/YvhphNF7brDGoX1FaDoOltsiCdcFgSRoSoIIqikL1ypbqgXWIiAB5duqgL2jVsqHF1pZNRkMH3h75n3uF5ZJuKrxvTv1F/XAwuGldYAxVkweFf1Usw0WvhXAM50EG9Lmo31BYDwEMulQnHJ2FEiApQeOwYSZPeIe/ffwEwhoQQNH4cXr16ValLMil5KXxz4Bt+Pvpz0boxET4RRevGGPVyr0GlMuerPUD2L4Sjq8BaeOG1Om3VEZDIu8Gnai8VIGoeCSNC2Fn2mjXEPTcKLBZ1Qbthw9QF7dyq1jX65SeW8/a/bxeFkOZ+zRnWahi96vWSdWMqW34GbJ0N/06H/PQL2wOanusFMhD8q9ZomxAXkzAihB3l79lD/NjnwWLB44Ye1P7fa1VuQTuzzcwH2z5g3uF5ALQOaM2TbZ6ke2j3KjWqUy3kpqkBZOusCw3JfMLU8NFqEAS3lF4golqQMCKEnZhiYoh9aiRKQQGeN9xA3c+noXOqWv/FUvJSeOGfF9iZrK7I+mTrJ3mqzVMY9FWvAVuVlp0Emz6D7V+BOU/dFtgcerygNiWTvw9RzVSt75RCOCjL2bPEDn8Ca3o6rpGRhH70YZULIruTdzN27VhS8lPwNHryTrd3uCn8Jq3LqlkyYmDjJ7Dzuwv3g9RpCz1eVGfDSEMyUU1Vre+WQjggW0EBcU+NxHT6NMbQUMK+mIHew0PrskpMURTmH5nPe9vew2Kz0Mi3EVNvnEp9n/pal1ZzpJ2ADR/Bnp8urJAb1hF6vASNbpZLMaLak5itEYvFwquvvkpERARubm40aNCAt956C5vNpnVpohQUq5WEF18if/du9D4+hM2ehVNgoNZllViBpYBXN77KpC2TsNgs9Knfhx9u+0GCSGVJPgSLhsG09rDrezWIRNwAD/8Kj62Exr0kiIgaQUZGNPLee+/xxRdf8M033xAZGcn27dt59NFH8fHxYdSoUVqXJ0rozHvvkb16NTpnZ8I+n4ZLgwZal1RicdlxjF07lkPph9Dr9IyNGstDLR6Sm1QrQ8JuWP+BulLueY37qPeEhHXQrCwhtCJhRCObN2+mf//+3H777QDUr1+fH3/8ke3bt2tcmSiptK+/5uy33wEQ8t67uLdvr3FFJbcpfhMvrX+JzMJMarnU4v0b3qdjnY5al1X9xWyBde/D8dUXtjXvp4aQOm20q0sIjZXpMs306dOJiIjA1dWVqKgo1q9ff9X9f/jhB9q0aYO7uzt16tTh0UcfJS0trUwFX4uiKNjy8jT5pShKievs1q0bf/31F0ePHgVgz549bNiwgdtuu61C/lyEfWX9sZLk96YAEPTii3j37atxRSWjKAqz985mxJ8jyCzMpKV/S36+82cJIhVJUSD6H/j6DviqtxpEdHpodQ+M/BeGfCdBRNR4pR4ZmT9/PqNHj2b69Ol07dqVmTNn0rdvXw4ePEh4ePgl+2/YsIGHHnqIqVOncueddxIfH8+IESMYNmwYS5YsscsXcTElP58j7aLsftySaLpzBzr3kq3N8fLLL5OZmUmzZs0wGAxYrVYmTZrEfffdV8FVivLK27mThJdeAkWh1v334/fYo1qXVCI5phwmbJjAmtg1AAxsPJDxHcdLK/eKoihwbLU6EhK3Vd2mN0Lb+6DraGlSJsRFSh1GPvroIx5//HGGDRsGwMcff8zKlSuZMWMGkydPvmT/f//9l/r16/Pcc88BEBERwZNPPsmUKVPKWXrVNn/+fL7//nvmzZtHZGQku3fvZvTo0YSEhPDwww9rXZ64gsLok8Q9NRLFZMLz5psJfmV8lbjH4kTGCUb/PZpTWacw6o1M6DiBgU0Gal1W9WSzqWvGrHsfkvaq2wwuEPUwdHkOfMO0rU8IB1SqMGIymdixYwfjxo0rtr13795s2rTpsu/p0qULEyZMYMWKFfTt25fk5GQWLlxYdK/E5RQWFlJYeGHNhaysrBLXqHNzo+nOHSXe3550pWj3/eKLLzJu3DjuvfdeAFq1asXp06eZPHmyhBEHZUlNJfaJJ7BmZuLapjWhH7yPzuD4zadWn17NqxteJc+SR7B7MFNvnEqrwFZal1X9WC3qyrnrP4CUw+o2owdc/xh0fga8amtbnxAOrFRhJDU1FavVSnBwcLHtwcHBJCUlXfY9Xbp04YcffmDIkCEUFBRgsVjo168fn3322RXPM3nyZN58883SlFZEp9OV+FKJlvLy8tD/p4GRwWCQqb0OypaXR+xTIzHHxWEMDydsxgyHX2vGYrPw6a5Pmbt/LgAdandgSo8p+LvJKq52ZTHB3vlqn5D0aHWbizd0fBI6PiWr5gpRAmW6gfW/w9KKolxxqPrgwYM899xzvPbaa+zYsYM//viDkydPMmLEiCsef/z48WRmZhb9io2NLUuZDu3OO+9k0qRJ/Pbbb5w6dYolS5bw0Ucfcdddd2ldmvgPxWIh/vkXKNi3D4OvL+GzZuLk56d1WVd1tuAsI/4cURREHol8hJm3zJQgYk/mAnXxus/awbJn1CDi5gc9X4XR+9TfJYgIUSKlGhkJCAjAYDBcMgqSnJx8yWjJeZMnT6Zr1668+OKLALRu3RoPDw+6d+/OxIkTqVOnziXvcXFxwcWlet9U99lnn/G///2PkSNHkpycTEhICE8++SSvvfaa1qWJiyiKQtKkSeT8/Tc6FxfqzpiOc/36Wpd1VQdSDzBm7RgScxNxc3Ljra5vcWv9W7Uuq/oozIEdc9W1Y3LOqNs8g6HLsxD1KLh4alufEFVQqcKIs7MzUVFRrF69uthP8KtXr6Z///6XfU9eXh5O/1mjw3DuOntppsJWN15eXnz88cd8/PHHWpciriJ9zhwyfvwJdDpCPngf9+uu07qkq1pybAkT/52IyWainnc9Pr7xYxrVaqR1WdVDQaa6eu7m6ZCfrm7zrgvdRsN1D4DRsS/bCeHISj2bZuzYsTz44IO0b9+ezp07M2vWLGJiYoouu4wfP574+Hi+/fZbQL0cMXz4cGbMmEGfPn1ITExk9OjRdOjQgZCQEPt+NULYUeavv5H8wYcABI8fj/ctt2hc0ZWZrCbe3fouC44uAODGsBt5p9s7eDl7aVxZNZCbBltmwJZZUJipbqsVAd2fh9ZDwMlZ2/qEqAZKHUaGDBlCWloab731FomJibRs2ZIVK1ZQr149ABITE4mJiSna/5FHHiE7O5tp06bx/PPP4+vrS8+ePXnvvffs91UIYWe5W7eSOH48AH4PP4zfQw9qXNGVJeUm8fza59mbuhcdOp5u+zTDWw9Hr5Olp8ol+wxs/gy2fQXmXHVbYDPo/gJE3gUGaWAthL3olCpwrSQrKwsfHx8yMzPx9vYu9lpBQQEnT54s6ghbndWkr1VLhcePc2ro/diysvDq04fQqR+hc9Cl27clbeOFf14gvSAdb2dv3u3+Lt3rdte6rKotIxY2fgI7vwXruRYDtVtDjxeh2R3goP8WhHBEV/v8vphEeyEuYk5OJuaJJ7BlZeHWrh0hU95zyCCiKArfHfyOj3Z8hFWx0rRWU6beNJUwL2moVWZpJ2DDVNjzo7p6LkDd66HHS9D4Flk9V4gKJGFEiHOsObnEjhiBJSER5/r1qfv5NPQOOKsrz5zHG5ve4PdTvwNwR4M7eK3za7g5yQ2UZZJ8GNZ/CPsXgnKuz0/97upISEQPCSFCVIJqE0aqwNWmcqsJX6NWFLOZ+DFjKDx4CIO/P2GzZ+FUq5bWZV0iJiuGUX+P4njGcZx0Trxw/QsMbTa0SrSkdziJe2DdB3BoOXDu/1ajW9QVdMM7aVqaEDVNlQ8jRqMRUKcQuzl4R8zyysvLAy58zcI+FEUh6a23yF2/Hp2bG2FfzMA5zPEud6yLW8e4dePINmcT4BbAhzd8SLvgdlqXVfXEblVDyLGVF7Y1v1OdHRPi2FO3haiuqnwYMRgM+Pr6kpycDIC7u3u1+ylRURTy8vJITk7G19e3qE+LsI+0L74gY8FC0OsJ/fBD3Fo51rotNsXGzD0zmb5nOgBtA9vy4Y0fEuQepHFlVYiiwKkN6uJ1J/9Rt+n00HIgdBsLwS20rU+IGq7KhxGA2rXVBajOB5LqytfXt+hrFfaR8csvpHzyKQC1//cqXj1v0rii4jILM3llwyusi1sHwL1N7+Wl61/CaJDRsRJLPgy/jYXTG9Xneidoc68aQvwbalubEAKoJmFEp9NRp04dgoKCMJvNWpdTIYxGo4yI2Fnupk0kvvo/APyHD6PWffdpXFFxR88eZfTfo4nNjsXF4ML/Ov2P/o0u3+lYXIbVrE7R/ec9sJrA4ALtHoSuo8A3XOvqhBAXqRZh5DyDwSAf2KJECo4cIe7Z58Biwfv22wkcM0brkopZEb2CNza/Qb4ln1DPUKbeOJXm/s21LqvqSNwLS0dC0j71eePecMdU8KmrbV1CiMuqVmFEiJIwJyUR+8ST2HJzcb/+eupMfsdheomYbWam7pjKdwe/A6BLSBfe6/4evq6+2hZWVVgK1ftCNkxVe4W4+kLf99S27dXsXjIhqhMJI6JGsWZnE/vEk1jOnMG5UUPqTvsMvbNjrC2Smp/KC/+8wI4zOwAY3mo4T7d9GoNeRvtKJG47LH0aUg6rz5v3g9s+AK/LryguhHAcEkZEjaGYTMQ99xyFR49iCAwgfOZMDD4+WpcFwJ6UPYz9eyzJ+cl4GD2Y1G0SN4ffrHVZVYMpD/6eBP9OV5uWeQSqISRygNaVCSFKSMKIqBEURSHxf6+Rt/lfdO7uhH3xBcbQUK3LQlEUFhxdwOStk7HYLDTwacDHN31MhE+E1qVVDac2wrJnID1afd56CNz6Lrj7aVuXEKJUJIyIGiH1s8/IXLoUDAbqfvIxbpGRWpdEobWQSf9OYsnxJQDcUu8W3u76Nh5GD40rqwIKs+HPN2Dbl+pzrxD1BtWmt2palhCibCSMiGrv7IIFpE6fAUCdN9/As7v2q9om5CQwZu0YDqYdRK/TM6rdKB6NfLTaNeyrEMf/guWjIDNWfd7uYej9Nrg6xiU3IUTpSRgR1VrO+vUkvfEmAAEjn8J30CCNK4LNCZt5ad1LZBRm4Oviy/s3vE+nOrIWyjXlZ8CqCbDre/W5bzj0+wwa3KhlVUIIO5AwIqqtgoMHiR81GqxWfPr3J+DZZzWtR1EU5h6Yyyc7P8Gm2Gjh34KpN04lxDNE07qqhMMr4NcxkJME6KDjk9Dzf+DiqXVlQgg7kDAiqiVzfDwxTz6JLS8Pjy6dqfP2W5peAsk15/K/jf9j9enVAAxoNIBXO72Ki8FFs5qqhNxU+P1l2L9Qfe7fCPpNg3qdta1LCGFXEkZEtWPNzCTmiSexpqTi0qQJoZ98gk7DXiLRmdGM+XsM0ZnROOmdGN9hPIObDJb7Q65GUeDAYljxIuSlqYvadXkWbhwPxuq9OrcQNZGEEVGt2Ewm4p55FtOJEzgFBxM2ayYGLy/N6vnr9F9M2DiBXHMuQe5BfHTjR7QJbKNZPVVCdhL89jwc/lV9HtQC+k+D0Cht6xJCVBgJI6LaUGw2Ese/Qt62beg9PQmbNROjRqscW21Wpu2expf71Kmn7YPb8/4N7xPgFqBJPVWCosDuebByPBRkqqvrdn8Buj8PTo7RJVcIUTEkjIhqI2XqVLJ++w2cnKj76Se4Nm2qSR255lzGrh3LpoRNADzU4iFGR43GqDdqUk+VkBGrTtc98Zf6vE5b6P851G6paVlCVHeZ+WZW7k9i2Z4Epg5pS6CXNvexSRgR1cLZH38kbbY6ClHn7bfx6NJFkzqsNisvr3uZTQmbcHNy480ub9I3oq8mtVQJNhvsmAurXwNTDhhc4Kbx0PlZMMi3JyEqQoHZyl+Hklm6O561R1IwWW0A/LY3gUe6atP9Wf63iyove83fJL09EYCA557F964BmtXyyc5P+CfuH1wMLnzZ+0taB7bWrBaHl3YClj0Hpzeoz8M6qqMhAY21rUuIashitbHxRBpLd8ez6sAZcgotRa81DfaiX9sQbonU5rI2SBgRVVz+vn3EP/882Gz4Dh5EwFNPaVbLL8d/Ye6BuQC83fVtCSJXYrPCli/gr7fBkg9Gd7j5degwHGSFYiHsRlEUdsZksGx3PL/uTSQt11T0WqivG/3bhtCvbQjNantrWKVKwoioskyxscSOeAolPx+P7t2p/dprmk2X3ZW8i7c2vwXAk62flEszV5JyBJY+DXHb1OcRPeDOT8FPFgYUwl6OJGWzdHc8y/YkEHc2v2i7n4czd7SuQ/+2IbQLr+VQ7QUkjIgqyXL2LLHDn8CaloZLi+aETp2KzqjNDaIJOQmM/ns0ZpuZW+rdwsi2IzWpw6FZzbDxE/jnPbCawNkL+kxU15VxoG+IQlRVsel5LN+bwLLdCRxOyi7a7uFsoE9kbfq1DaFrowCMBr2GVV6ZhBFR5dgKCoh7+hlMp07hFFKHsC++wOCpzUq3eeY8nl3zLOkF6TT3a87ErhPR6xzzP7tmEveqoyFJe9XnjXurK+z61NW2LiGquNScQlbsS2Tp7gR2nD5btN3ZoOfGpoH0axvCzc2CcXN2/MufEkZElaLYbCS8PI78nTvRe3kRPmsWxqAgTWqxKTbGrR/H0bNH8Xf159Oen+JudNekFodkKYR1H8CGj8BmAVdf6PsetB4ioyFClFFOoYVVB5JYujuBDcdTsdoUQP0v1bmBP/3bhnBrZB183KtWKwEJI6JKSX5vCtkrV6IzGqk7bRoujRppVstnuz7j79i/cdY782nPT6ntod2d6A4nbrs6GpJyWH3evB/c9gF4BWtblxBVUKHFytojKSzbk8CfB89QaLEVvda6rg/92oRwZ5sQgr1dNayyfCSMiCoj/dtvSf/mGwDqTJ6MR8cOmtWy/MTyou6qb3Z9U2bOnGfKg7XvwObPQbGBR6AaQiIHaF2ZEFWK1aawJTqNpbsT+H1/IlkFF6biNgjwoH/bUPq1DSEiQJtL1PYmYURUCVmrVnFm8rsABD4/Fp87bteslj0pe3hj0xsADGs1jDsa3KFZLQ7l1EZY9gykR6vPWw+BW98Fdz9t6xKiilAUhb1xmSzbk8DyPQkkZxcWvVbb25U729Shf9tQIkO8HWomjD1IGBEOL2/nLhJefAkUBd/77sV/2DDNaknKTWLUmlGYbCZuCruJZ697VrNaHEZhNvz5JmybrT73ClFvUG16q7Z1CVFFnEjJYenuBJbtjudUWl7Rdh83I7e1Uqfidqjvh15fvQLIxSSMCIdWePIkcSNHohQW4nnjjdSeMEGznwjOz5xJK0ijSa0mvNv9XZk5c2INLBsFmTHq83YPQ++3wdVH27qEcHCJmfn8uieRpXvi2R+fVbTd1ajnlha16d8mhB5NAnF2qhnfYySMCIdlSUsj9oknsWZk4NqyJaEffYjOSZt/sjbFxqsbX+Vw+mH8XP34rOdnNXvmTH4GrJoAu75Xn/uGQ7/PoMGNWlYlhEPLyDOxYl8SS3fHs/VUOoo6EQaDXkePxgH0bxvKLS2C8XCpeR/NNe8rFlWCLT+f2KdGYo6NxVi3LmFfzEDvrt2H//Td01l9ejVGvZGPb/qYEM8QzWrR3OEV8OsYyEkCdNDxSej5P3Dx1LoyIRxOnsnCn4eSWbY7nn+OpmC2KkWvdajvR7+2IdzWqg5+Hs4aVqk9CSPC4ShWK/EvvEjB3r0YfHwImzULp4AAzer5/eTvzNw7E4DXO7/OdUHXaVaLpnLT4PeXYP9C9bl/I+g3Dep11rYuIRyM2Wpj/bEUlu1OYNXBM+SZrEWvNa/jTf+26lTcUF83Dat0LBJGhENRFIUzk94h56+/0Dk7U3fGdFwaaLduyf7U/fxv4/8AeDTyUfo36q9ZLZpRFDiwBFa8CHmpoNNDl2fhxvFglG+mQgDYbArbT59l6e54VuxL5Gyeuei1cD93+rVRF6VrEuylYZWOS8KIcCjpX83l7Lx5oNMRMmUK7u3aaVbLmdwzPLfmOQqthdxQ9wZGtRulWS2ayU6C356Hw7+qz4NaQP9pEBqlbV1COABFUTiUmM3SPfEs351AQmZB0WsBns7c0TqE/m1DaBvmW+2m4tqbhBHhMHLWryf5/fcBCHrpJbxv7aNZLfmWfEb9PYqU/BQa+Tbi3e7vYqhpy9tHr4WfH4aCDNA7QfcXoPvz4FSzr20LkVVg5pdd8czbElNsUTpPFydubVmb/m1D6NzAHycHXZTOEUkYEQ5BsVpJnqIGEd8hQ/B75GHtalEUXtv4GgfSDuDr4sunPT/F07mG3Zy581v1JlWbBeq0gf7ToXZLrasSQjOKorA7NoN5W2JYvjeBArPakt3ZSU/PpkH0bxvCTc2CcDXWsB9a7ETCiHAIWStWUHjsGHpvb4LGjtF0SPOLvV/wx6k/cNI7MfXGqYR5hWlWS6Wz2eCvN2DjJ+rzVoPVm1SNVXfNCyHKI6vAzNJd8fzwn1GQxkGeDO0Yzt3X1a1yi9I5IgkjQnOK2UzKp58B4P/44xh8tGuYterUKqbvng7A/zr9j/a122tWS6Uz5cGSJ+DQcvX5DePgxnGywq6ocRRFYU9cJvO2nGb5nkTyzepsGGcnPXe0qsPQjuFE1asl94HYkYQRobmMRYsxx8Zi8PfH78EHNKvjYNpBJmyYAMADzR/g7sZ3a1ZLpctOgh/vhYRdYHCG/p9D63u0rkqISpVdYOaX3Qn8uCWGg4kXuqI2CvJkaIdw7m4Xiq+73DNVESSMCE3ZCgpIna6ORASMGKFZY7OUvBSeXfMsBdYCuoZ25fn2z2tShybOHIAf7oGsOHDzg3t/gHpdtK5KiEpxfnG6eVtiWLYnodgoyO3nRkHayyhIhZMwIjR19od5WJKTcQqpg+8QbX4SL7AUMOrvUSTnJRPhE8H7Pd7HSV9D/mscWw0LHgVTttrEbOjP4N9Q66qEqHDZBWaW7k7gx60xHEi4MArSMNCDoR3rMVBGQSpVDfmOKxyRNSeHtFmzAAh8+hn0zpX/H19RFF7f9Dr7Uvfh4+LDtJ7T8HKuIU2Jts5WO6oqNqjfHe75Ftz9tK5KiAq1Ny6jaBTkfGdUZyc9t7WszdCO9bi+voyCaEHCiNBM+tyvsWZm4tygAT79+2lSw5f7vmTFyRU46Zz46IaPCPcO16SOSmWzwsoJsGWG+rztA3DHVOkfIqqtnEILS3fH8+PWmGIr5DYI9GBoh3AGtqtLrRq+NozWJIwITVjOniV97lwAAp97TpPVeP86/Ref7voUgPEdx9OhTodKr6HSFebAosfh6B/q85tfh25jZMaMqJb2xWUyb+tplu6+aBTEoKdvq9oM7RBOhwg/GQVxEBJGhCbSZs3GlpeHa4sWePW+pdLPfyT9COM3jAfgvmb3cU/TGjBzJDMefhwCSfvAyRXu+gIi79K6KiHsKqfQwrJz94Lsi88s2n5+FOTudnVr/Aq5jkjCiKh05qQkzv7wAwCBY0aj01duy+TU/FSeWfMM+ZZ8OtXpxEvXv1Sp59dEwi6Ydy/kJIFHINz3E9StQT1URLW3Pz6TH7bEsGx3PLkXjYLc2rI2QzuG01FGQRyahBFR6VKnz0AxmXBrH4VHt26Vem6T1cTov0eTlJtEfe/6fHDDB9V/5szh32DRMDDnQWBzGDofatXTuiohyi230MKyPeooyN64i0ZBAjy4r0M4A6NkFKSqqObfhYWjMZ0+TcaiRQAEjR5dqT+pKIrCm5vfZE/KHrycvfis52f4uGjX7bXCKQps/hxWvQoo0LAnDP4aXKvx1yxqhP3xmczbGsPSXRdGQYwGHbe2rMPQDuF0aiCjIFWNhBFRqVI+mwZWKx49uuPevnIvE3x94GuWnViGQWfggxs+oL5P/Uo9f6WymmHFi7BDvUmY9o9B3/fBIP/lRdWUW2hh+blRkD0XjYJEBHhwX4cwBrari7+ni4YVivKQ70yi0hQcOULWb78B6qhIZVobu5apO6YC8NL1L9ElpBp3GC3IhAWPwIk1gA76TIJOI2XGjKiSDiSo3VGX7k4gp9ACqKMgfSLVe0E6N/CXUZBqQMKIqDQpn3wKioJX31txbdGi0s579OxRXl73MgoK9zS5h/ua3Vdp5650Z0/DvCGQcgiM7jBwDjS7TeuqhCiVPJM6CjJvayx7YjOKttf3d+e+DuEMipJRkOpGwoioFHm7dpGzZg3o9QQ++1ylnTe9IJ3n1jxHniWPDrU7MK7juOr7U1TsNvjpPshNAa866oyZkLZaVyVEiR1MyOLHrTH8siue7ItGQXpH1ub+DuF0auCPXl9N///WcGWaUzl9+nQiIiJwdXUlKiqK9evXX3X/wsJCJkyYQL169XBxcaFhw4Z89dVXZSpYVD2KopAy9WMAfO4agEuDiEo5r9lqZszfY4jPiSfMK4wPb/gQo95YKeeudPsXwzd3qEGkdisY9pcEEVEl5Jks/LwtlgGfb+S2T9fz3b+nyS60UM/fnXF9m7F5/M18PrQdXRoFSBCpxko9MjJ//nxGjx7N9OnT6dq1KzNnzqRv374cPHiQ8PDLt9K+5557OHPmDHPmzKFRo0YkJydjsVjKXbyoGvI2byZv61Z0RiOBI0dWyjkVReHtf99mZ/JOPI2eTOs5DV9X30o5d6VSFFj/Iax5W33epC8M/BJcPLWtS4hrOJSYxbwtxUdBnPTF7wWR8FFz6BRFUUrzho4dO9KuXTtmzJhRtK158+YMGDCAyZMnX7L/H3/8wb333kt0dDR+fmVbhCsrKwsfHx8yMzPx9vYu0zGENhRF4dQ9QyjYt49aDz5I7QmvVMp5vz3wLe9vfx+9Ts/nN39Ot9DK7WdSKSwm+HU07FYbyNFpJPSeCHqDpmUJcSVWm8Lqg2f4cn0020+fLdoe7nfhXpBAL7kXpDop6ed3qUZGTCYTO3bsYNy4ccW29+7dm02bNl32PcuWLaN9+/ZMmTKF7777Dg8PD/r168fbb7+Nm5tbaU4vqqCcv/6iYN8+dO7uBDz5RKWcc33cej7c8SEAL7R/oXoGkbx0mP8gnN4AOj30nQIdhmtdlRCXlW+ysnBHLHM2nORUWh6gjoL0jgxmaId6dGkooyA1XanCSGpqKlarleDg4GLbg4ODSUpKuux7oqOj2bBhA66urixZsoTU1FRGjhxJenr6Fe8bKSwspLCwsOh5VlbWZfcTjk2xWkn55BMA/B56EKeAgAo/54mME7y07iVsio2BjQfyQPMHKvyclS7tBMy7B9KOg7OX2siscS+tqxLiEinZhXy3+RTf/Xuas3lmAHzcjDzQKZyHO9cnyNtV4wqFoyjTbJr/zkZQFOWKMxRsNhs6nY4ffvgBHx+18+NHH33EoEGD+Pzzzy87OjJ58mTefPPNspQmHEjWr79SeOw4em9v/B97rMLPl1GQwTN/PUOOOYeo4CgmdJxQ/WbOnN4EPw2F/LPgE6a2dg+O1LoqIYo5npzDnA3RLNoZj8liAyDMz43Hu0YwuH0YHi4ykVMUV6p/EQEBARgMhktGQZKTky8ZLTmvTp06hIaGFgURUO8xURSFuLg4GjdufMl7xo8fz9ixY4ueZ2VlERYWVppShcYUk0nttgr4DxuGoYLv9TFbzYz9ZyxxOXGEeoYy9capGA3VbObMnvmw7BmwmiCknTp11+vy/++EqGyKorDlZDqz10Xz1+Hkou1tw3x5okcD+kTWxiCXYsQVlCqMODs7ExUVxerVq7nrrgtLj69evZr+/ftf9j1du3ZlwYIF5OTk4Omp3uF/9OhR9Ho9devWvex7XFxccHGRm5iqsoxFizDHxWEICMDvgfsr9FyKovDO1nfYlrQNdyd3pvWcRi3XWhV6zkqlKLB2Mvzznvq8eT+4ayY4u2tblxCAxWrj9/1JzF4fXbRYnU4HtzQP5okeDYiqV6v6jVAKuyv1WNnYsWN58MEHad++PZ07d2bWrFnExMQwYsQIQB3ViI+P59tvvwVg6NChvP322zz66KO8+eabpKam8uKLL/LYY4/JDazVlC0/n9Tp6myrgBEj0LtX7IfmvMPzWHh0ITp0TOkxhUa1GlXo+SqVuQCWPg37F6rPu42Bnq+BvkwtgoSwm5xCC/O3xfLVhpPEZ+QD4OKkZ1BUXR7vFkGDQJleLkqu1GFkyJAhpKWl8dZbb5GYmEjLli1ZsWIF9eqpS5InJiYSExNTtL+npyerV6/m2WefpX379vj7+3PPPfcwceJE+30VwqGcnTcPS0oKxtBQat0zuELPtSl+E1O2TQFgbNRYbgi7oULPV6lyU9X7Q2K3gN4J7pgK7R7SuipRwyVlFjB300nmbYkhu0DtD+Lv4cxDnevzQKdwadMuyqTUfUa0IH1Gqg5rdjYnet2CNTOTOpMn43vXgAo7V3RmNA/89gDZ5mz6NezHxK4Tq89wcMoR+GEwZJwGVx+45ztoUI2ClqhyDiVmMXt9NMt2J2CxqR8bDQI8GNa9AXe3C8XVKP1txKUqpM+IENeSPncu1sxMnBs2xKffnRV2nszCTJ7961myzdlcF3Qdr3d+vfoEkei1MP8hKMyEWvVh6AIIbKJ1VaIGUhSF9cdSmb0+mvXHUou2d4jw44nuDejZLEj6gwi7kDAi7MaSlkba198AEPjcc+gMFfOTktlm5vl/nicmO4Y6HnWYeuNUnA3OFXKuSrfzW/h1DNgsENYJ7p0HHv5aVyVqGJPFxvI9CcxeH83hpGwA9Dq4rVUdhndvQJswX20LFNWOhBFhN2mzZqPk5eEaGYlX71sq7DxTtk5hS+IW3Jzc+KznZ/i7VYMPa5sN/noTNn6sPm81GPpNA6M0hRKVJzPfzLwtMXy96SRnstTGk+7OBoZcH8ZjXSMI85MZXKJiSBgRdmFOTOTsjz8CEDh6dIVdMpl/eD4/HfkJHTre7f4uTf2aVsh5KpUpD5Y8AYeWq89vGAc3jlPnRwpRCWLT85i78RTzt8WQa7ICEOTlwqNdIxjaIRwf92rWs0c4HAkjwi5Sp89AMZlwv/56PLp1rZBz/Jv4L5O3qosxPtfuOXqG96yQ81Sq7CT48V5I2AUGZ+j/ObS+R+uqRA2xNy6DWeuiWbEvkXP3pNI02IvhPRrQr00Izk4yhVxUDgkjotxMp06RsXgxAIFjKmZU5HTWaZ5f+zxWxcodDe7g8ZaP2/0cle7MAfjhHsiKAzc/uPcHqNdF66pENWezKaw5nMzs9dFsOZletL174wCGdW9Aj8YB1edmcFFlSBgR5Zby6WdgteJ5ww24t2tn9+NnmbJ45q9nyDJl0TqgNW90eaPqf7M8thoWPAqmbPBvBEN/Bv+GWlclqrECs5Ulu+KZvT6a6JRcQF05t1+bEIZ1b0CLEGmbILQjYUSUS8Hhw2StWAFA4OhRdj++xWbhxX9e5FTWKYLdg/mk5ye4GKp4U6Wts+H3l0CxQf3uMOQ7cKtG7euFQ0nPNfH9v6f5dvMpUnNMAHi5ODG0UziPdKlPHR/phC20J2FElEvKx58A4H1bX1ybN7f78T/c/iGbEjYVzZwJcAuw+zkqjc0KKyfAFrVVPm0fULuqOlWTacnCoZxMzWXOhmgW7oijwKyunBvq68ajXesz5PowvFzlplThOCSMiDLL27mLnLVrwWAg4Nln7X78hUcX8v2h7wGY1G0Szf3tH3YqTWEOLHocjv6hPr/5dXWdmap+uUk4FEVR2HH6LLPWRbP60BnO99duFerD8B4NuK1lbZwMclOqcDwSRkSZKIpCytSpAPjefRcuERF2Pf62pG1M+ncSAE+3fZpb6lVc35IKlxkPPw6BpH3g5Ap3fQGRd137fUKUkNWmsOpAErPWR7MrJqNo+83NghjeowEdI/yq/n1WolqTMCLKJHfTJvK2bUNnNBIwcqRdjx2bHcuYtWOwKBb61u/Lk62ftOvxK1XCLvjxPshOBI9AuO8nqNte66pENZFnsrBgexxzNpwkJj0PAGeDnrvbhTKsewSNgrw0rlCIkpEwIkpNHRX5GIBaQ+/DWKeO3Y6dY8rh2b+eJbMwk0j/SN7q+lbV/Ynu8G+waBiY8yCwOQydD7XqaV2VqAaSswv4dtNpvvv3NJn5ZgB83Y081KkeD3auT6BXFb/JW9Q4EkZEqWWvXk3B/v3o3N3xf+IJux3XarPy0rqXOJF5giC3ID7t+SmuTlWwHbqiwObPYdWrgAINe8Lgr9XVd4Uoh+PJ2cxaF80vuxIwWdWbUuv5uzOsWwQDo+ri7izf0kXVJP9yRakoVispn3wKgN/DD+Hkb791YeYemMv6+PW4GFz4tOenBLkH2e3Yleqf92Ct2imW9o9D3ylgkP9qouyOncnmk7+O8du+xKKbUqPq1WJ49wbc0iIYg6ycK6o4+Q4pSiVz+XJMJ06g9/HB/7HH7HbchJwEZu6ZCcCEjhOIDIi027Er1aZpF4JIrzeh6yiZMSPK7HhyDp/+dYzlexOKQkjvFsE8eUMDour5aVucEHYkYUSUmGIykfrZNAAChg/D4GW/m+Pe2/oeBdYCooKjGNBogN2OW6l2fA2rJqiPe74K3UZrWY2owk6k5PDZX8dYtiehaM2YPpHBjLq5iXRKFdWShBFRYmcXLsQcH48hMIBa999vt+Oui1vHmtg1GHQGJnScUDVvWN27AJaPVh93HQ3dX9CyGlFFRafk8Nma4yzdHV8UQnq3CGZUr8ZEhsg9R6L6kjAiSsSWl0fqDLVzaMBTT6F3s08L6QJLAZO3qJc1Hmj+AI1rNbbLcSvV4d9gyZOAAtcPg15vyKUZUSqnUnP5dM0xftl1IYT0ah7M6F6NaRkqIURUfxJGRImk//AD1pRUjKGh1Bo0yG7H/Wr/V8TlxBHkFsRTbZ+y23ErzYm/YcEjoFih9b3Q930JIqLETqfl8tma4yzZFY/1XAq5uVkQo3s1oVVdCSGi5pAwIq7JmpVF2pdzAAh49hl0zvZZSyU2K5Y5+9TjvtjhRTyMHnY5bqWJ2QI/DQWrCZrfCf0/B7202hbXFpOWx2drjrH4ohDSs1kQo25uTJswX22LE0IDEkbENaXNnYstMxPnRg3xufNOuxxTURTe2foOJpuJTnU60adeH7sct9Ik7IYfBqsNzRr1goFzZPquuKbY9DymrTnOop1xWM6FkBubBjK6VxPaSggRNZh89xRXZUlLI/2bbwEIHDUKncFgl+OuiVnDhvgNOOmdeKXjK1XrptXkw/D93VCYCeFd4J7vwEk6Xoori03PY/ra4yzYfiGE9GgSyOhejWkXXkvj6oTQnoQRcVVps2ah5OXh2rIlXr162eWYeeY83t32LgCPRj5KhI99F9mrUOkn4bsBkJcGIdepLd6d3bWuSjio+Ix8pq05zsIdsZitagjp3jiA0b2aEFVPQogQ50kYEVdkTkjg7LwfAQgcM9puoxcz984kKTeJEI8QhrcebpdjVoqsBPi2v7roXWBzeGAxuErPB3GphIx8Pv/7OD9vvxBCujUKYHSvxrSvL83KhPgvCSPiilKmT0cxm3Hv0AGPLl3scszojGi+PaBe9hnXYRxuTvaZIlzhclPVIJJxGvwawEO/gLt8qIjiEjPzmf73CeZviy1aO6ZrI39G92rC9RJChLgiCSPisgqjT5K55BfAfqMiiqIwacskLIqFG+rewE3hN5X7mJUiP0O9NJN6FLzrwkNLwau21lUJB5KUWcD0tcf5aeuFENK5gT+jezWmYwP7rd8kRHUlYURcVuq0z8BqxfOmm3C/7jq7HPP3k7+zNWkrLgYXxnUYZ5djVrjCHHXWTNI+8AhUg4hvuNZVCQdxJquAGWtPMG9rDCaLGkI6RvgxulcTOjeUECJESUkYEZcoOHSIrBW/AxA4epRdjpljyuGD7R8AMKzVMOp61bXLcSuUuUDtIxK3FVx94cFfIKCR1lUJB5CcVcCMf04wb0sMhedCSIf6foy+pTFdGgZoXJ0QVY+EEXGJlI8/AcD79ttxbdrULsf8fPfnpOSnEO4VzqMtH7XLMSuU1ax2Vj35Dzh7qjer1m6pdVVCY8nZBXyxNpoftpwuCiHt69VizC1N6NLQv2pNURfCgUgYEcXk7dxJzj//gMFA4LPP2OWYR9KP8ONhdVbOKx1fwcXg4D05bFZ1rZmjv4OTK9z3E9SN0roqoaGU7EJm/nOC77ecpsCshpCoerUY06sJXRtJCBGivCSMiCKKopD80UcA+N59N87165f7mDbFxsR/J2JVrNxS7xa6hnYt9zErlKLAr6Nh/yLQG9WGZhHdta5KaCQ1p5BZ66L5dvOpohByXbgvY3o1oXvjAAkhQtiJhBFRJHfDRvK370Dn7EzA0yPtcsylx5eyO2U3bk5uvHT9S3Y5ZoVRFFg5AXZ+Czo9DJwNTXprXZXQQFpRCDlNvtkKQNswX8bc0oQeEkKEsDsJIwJQR0VSpk4FoNbQoRhrl3/qamZhJlN3qMd8qs1T1PZw8Omwa9+Ffz9XH/ebBpF3aVuPqHTpuaaikZA8kxpC2tT1YfQtTbixSaCEECEqiIQRAUD2qtUUHDyI3t0d/yfs0xX1052fcrbwLA19GvJAiwfscswKs+kz+EdtUU/f9+G6+7WtR1Sqs7kmZq+P5ptNp8g9F0Ja1/VhdK/G3NQ0SEKIEBVMwohAsVpJ+USdQeP3yCM4+ZW/U+T+1P0sOLoAgAmdJmDUG8t9zAqz/StY9ar6+ObXoOMT2tYjKk1GnhpCvt54IYS0DPVmTK8m9GwmIUSIyiJhRJC5dBmm6GgMPj74PfpIuY9ntVmZ+O9EFBRub3A719e+vvxFVpS9P8OvY9XH3cZA9+e1rUdUisw8M19uiGbuxlPkFFoAiAzxZnSvJvRqLiFEiMomYaSGs5lMpE6bBoD/E8MxeHmV+5iLji3iQNoBPI2evND+hXIfr8Ic+hWWjAAUuH443Py61hWJCpaZb2bOhpPM3XCS7HMhpHkdb0b3akzvFsESQoTQiISRGi7j5wWYExJwCgyk1tCh5T5eekE6n+xUL/k8c90zBLg5aDfKE2tg4aOgWKHNUOg7BeSDqNrKzDfz1YaTfLXxJNkFaghpVtuL0b2a0LtFMHq9/N0LoSUJIzWYLS+P1C++ACDg6ZHo3cq/gu7UHVPJMmXRzK8ZQ5oOKffxKsTpzfDT/WA1QfN+0O8z0Ou1rkpUgEKLle82n+azNcfJzDcD0DTYi9G9GtMnsraEECEchISRGiz9+x+wpqZiDAvD9+67y328Xcm7+OX4LwBM6DgBJ70D/vNK2AXz7gFzHjTqBQPngMEB6xTloigKy/cm8v7Kw8Sm5wPQOMiT0b2a0LelhBAhHI18F66hrFlZpH35JQCBzz6Dztm5XMez2CxM/HciAHc3vpu2QW3LW6L9JR+G7+6Gwiyo11XtrupUvq9bOJ4t0Wm8s+IQe+IyAQj2duH5W5oyMKouBgkhQjgkCSM1VNqcr7BlZeHSuBHet99e7uP9ePhHjp49io+LD6PbjS5/gfaWHg3f9of8dAhpp6434+yudVXCjk6k5PDu74dZffAMAB7OBkbc0JDHu0fg7izf6oRwZPI/tAaypKaS/u23AASOGoXOYCjX8ZLzkvl8t9q5dFS7UdRyrVXuGu0qMx6+6Q85SRDUAh5YBK7eWlcl7CQ1p5CP/zzKj1tjsdoUDHod914fxuheTQj0cvBFGYUQgISRGil15iyU/HxcW7fG8+aby328D7Z/QK45l1YBrRjYeKAdKrSjnBR1RCQzBvwawoO/gHv5m7oJ7eWbrMzZEM2MtSeKGpb1ah7MuL7NaBTkqXF1QojSkDBSw5jj48n46ScAgsaMLndfhS2JW/j95O/o0DGh0wT0OgealZJ/Fr67C9KOgXddeGgpeAVrXZUoJ6tNYdHOOD5adZSkrAJAXT9m/G3N6dTAX+PqhBBlIWGkhkmZPh3FbMa9Uyc8Oncu17HMVjOTtkwC4J6m9xDpH2mPEu2jMAd+GAxn9oFHEDy8DHzDtK5KlNO6oym8s+IQh5OyAQj1deOlW5tyZ+sQmSEjRBUmYaQGKYw+SeaSXwAIGj2q3Mf79uC3nMw8iZ+rH89e92y5j2c35gL46T6I2wZuteChX8C/odZViXI4lJjFOysOsf5YKgDerk4807MRD3Wuj6uxfPc8CSG0J2GkBkn59FOw2fDs2RO3tm3LdazEnERm7p0JwNiosfi4+NihQjuwmmHBw3ByHTh7qjerBjvQiI0olaTMAj5cdYSFO+NQFDAadDzUuT7P3NSIWh4yLVuI6kLCSA2Rf+AA2X/8ATodgaPKPyoyZdsU8i35tAtqR7+G/exQoR3YrLD4CTj6Bzi5wtD5EBqldVWiDHIKLXyx9gRfboimwGwD4PbWdXipT1Pq+XtoXJ0Qwt4kjNQQKZ+o68V433EHrk2blOtY6+PW82fMnxh0BiZ0muAYi4vZbLB8FBxYDHojDPkB6nfTuipRSmarjZ+2xfLJn0dJzTEBcH39WrxyW3OuC3ewKeNCCLuRMFID5G3fTu669eDkROAzT5frWIXWQiZvnQzA/c3vp0mt8gUbu1AUWPkK7PoOdHoYNAca99K6KlEKiqKw+uAZ3v3jMNEpuQBEBHgwrm8zWU1XiBpAwkg1pygKyVM/BsB34ECc69Ur1/G+2vcVsdmxBLkFMbLtSDtUaAd/vwNbZqiP+38OLfprW48olT2xGUxacYitJ9MB8PNwZnSvxtzXIRyjwYGmigshKoyEkWoud8MG8nfsQOfsTMDIp8p1rNisWL7cp65n8+L1L+JhdIBr9xs/gXVT1Me3fQBth2pbjyix2PQ8pqw8wvI9CQC4OOl5vFsEI25siLerUePqhBCVScJINabYbCRPnQpArfvvxxhc9oZfiqIweetkTDYTnep0ok/9PvYqs+y2zYHVr6mPb34dOgzXth5RIpl5Zqb9fYxvNp3GZLWh08Hd19Xl+d5NCPF107o8IYQGJIxUY9mrVlF48BB6Dw/8nyjfB/Wa2DWsj1+Pk96JVzq+ov01/D3z4bfn1cfdn4fuY7WtR1xTocXKd5tP89ma42TmmwHo1iiA8bc1IzLEQaaGCyE0IWGkmlIsFlI++RQAv0cfxalW2Wci5JnzeG/rewA8EvkIET4RdqmxzA4th1+eAhTo8CT0/J+29YirUhSFX/cmMmXlYWLT8wFoGuzF+NuacUOTQO2DrRBCc2W6O2z69OlERETg6upKVFQU69evL9H7Nm7ciJOTE23L2XBLXFvm0mWYTp7E4OuL3yMPl+tYs/fNJjE3kRCPEJ5o/YSdKiyj43/CgkdBsULb++HWd0E+zBzW1pPpDJi+iWd/3EVsej5BXi68N7AVK0Z158amQRJEhBBAGUZG5s+fz+jRo5k+fTpdu3Zl5syZ9O3bl4MHDxIeHn7F92VmZvLQQw9x8803c+bMmXIVLa7OZjKR8vk0APyfeAKDZ9lXMI3OjObrA18D8HKHl3Fz0vCa/ulN8NMDYDNDiwHQ7zPQy2wLR3QiJYd3fz/M6oPq/3V3ZwMjbmjIsO4RuDvLgKwQojidoihKad7QsWNH2rVrx4wZM4q2NW/enAEDBjB58uQrvu/ee++lcePGGAwGfvnlF3bv3l3ic2ZlZeHj40NmZibe3t6lKbdGSv/ue85MmoRTUBANV61E7+papuMoisLw1cPZkriFHnV7MK3nNO1+ko3fCd/0A1M2NO6tNjVzknbgjiY1p5BP/jzGvK0xWG0Keh3c2yGc0b0aE+RVtn+HQoiqq6Sf36X6EcVkMrFjxw7GjRtXbHvv3r3ZtGnTFd83d+5cTpw4wffff8/EiROveZ7CwkIKCwuLnmdlZZWmzBrNlptL6hdfABAwcmSZgwjAH6f+YEviFlwMLozrME67IJJ8CL4fqAaRet3gnm8liDiYfJOVrzaeZMbaE+QUWgC4uVkQ4/o2o3Gwl8bVCSEcXanCSGpqKlarleD/TBENDg4mKSnpsu85duwY48aNY/369Tg5lex0kydP5s033yxNaeKc9O++x5qWhjE8HN+Bd5f5ODmmHN7f9j4Aj7d6nDCvMHuVWDppJ+Db/pCfrq4zM/QnMMr0T0dhtSks3hnHh6uOkpRVAECrUB9eua05nRv6a1ydEKKqKNPF2//+hKwoymV/arZarQwdOpQ333yTJk1K3jZ8/PjxjB17YapmVlYWYWEafRhWIdbMTNLmzAEg8Nln0RnL3jhq+p7ppOSnEOYVxmMtH7NXiaWTGQffDoCcMxDcEu5fCC7yU7ajWH8shXdWHOZQojpyGerrxku3NuXO1iHo9XJjqhCi5EoVRgICAjAYDJeMgiQnJ18yWgKQnZ3N9u3b2bVrF8888wwANpsNRVFwcnJi1apV9OzZ85L3ubi44OLiUprSBJA25yts2dm4NGmC9+23lfk4R88eZd6heQC80vEVXAwa/F3kJKsjIpkx4NcQHlwC7n6VX4e4xKHELCb/fph1R1MA8HJ14pmbGvFwl/q4Gg0aVyeEqIpKFUacnZ2Jiopi9erV3HXXXUXbV69eTf/+l64H4u3tzb59+4ptmz59OmvWrGHhwoVERGjcr6IasaSkkP7ddwAEjh6FroyzTBRFYdK/k7AqVnqF96JbqAYr3+afhe/ugrTj4BMGDy0Fz6DKr0MUk5RZwIerjrBwZxyKAkaDjgc71efZno2o5SH38Aghyq7Ul2nGjh3Lgw8+SPv27encuTOzZs0iJiaGESNGAOollvj4eL799lv0ej0tW7Ys9v6goCBcXV0v2S7KJ3XmLJT8fFzbtMbzppvKfJxlJ5axM3knbk5uvNzhZTtWWEKF2fD9IDizHzyD1SDiK5fotJRTaGHmPyeYvT6aArMNgNtb1eGlW5tSz98B1icSQlR5pQ4jQ4YMIS0tjbfeeovExERatmzJihUrqHduNdjExERiYmLsXqi4MlNcPGfnzwcgaMyYMs96ySzM5KMdHwEwos0IanvUtluNJWLOhx/vg/jt4FYLHvwF/BtWbg2iiNWm8OPWGD7+8yipOSYAourV4pXbmhNVr+wdfYUQ4r9K3WdEC9Jn5OqS3nqbs/Pm4d65E/Xmzi3zcSb+O5H5R+bTwKcBC+9ciNFQiSunKgosGgb7F4KzFzy8DELbVd75RTHbT6Xz2tIDHDx3c2p9f3fG9W1Gn8ja0jVVCFFiFdJnRDgeW34+mcuXA+A/bFiZj3Mg7QA/H/kZgFc7vVq5QQRg48dqENE7wX3zJIhoJDm7gHd/P8zinfEAeLs6MeaWJtzfsR7OTtLtVghRMSSMVHHZq1Zhy87GGBqKR+fOZTqG1WZl4uaJKCjcFnEb19e+3s5VXsPRlfDnub4yfd+DiB6Ve36B2Wrjm02n+PjPY+QUWtDpYEj7MF7s0xR/T5nZJoSoWBJGqriMBQsB8B00sMwzaBYdW8T+tP14Gj15of0L9izv2lKOqJdnUCDqUbi+7KM7omw2n0jj9WX7OXomB4DWdX14q39L2ob5aluYEKLGkDBShRVGnyRv+3bQ6/G5u2zdVtML0vlk5ycAPN32aQLdA+1Z4tXlZ6g3rBZmQXgX6Dul8s4tSMzM550Vh1m+JwGAWu5GXr61Gfe0D5OmZUKISiVhpArLWKSOinj26IHxMk3nSuLjHR+TZcqiaa2m3NvsXnuWd3U2Kyx6HNJPqL1EZL2ZSmOy2Jiz4SSfrTlGnsmKXgf3d6zH872b4OsufwdCiMonYaSKUkwmMn9ZCoDv4EFlOsbu5N0sOb4EUG9addJX4j+HP1+H43+CkxvcOw88K3FEpgZbdzSFN5YdIDo1F4B24b681b8lLUN9NK5MCFGTSRiporL/Xos1LQ1DYACePUp/w6fFZmHiv+oKygMaDaBtUFs7V3gVe36CTZ+pjwdMhzqtK+/cNVTc2Twm/nqIPw6oSzkEeLowvm8z7rouVC7JCCE0J2GkispYeO7G1bvuLtOCePOPzOfI2SN4O3szJmqMvcu7svgdsOw59XH3F6Bl2VcWFtdWYLYya10009cep8Bsw6DX8XDn+oy+pTHerpU8fVsIIa5AwkgVZE5IIHfDBgB8B5b+wzwlL4Vpu6YBMKrdKPxcK2kBuuwk+Ol+sBZCk75w04TKOW8N9dehM7y5/CAx6XkAdIzw483+kTSrLY0DhRCORcJIFZSxaDEoCu4dO+J8rg1/aXy440NyzDm09G/JwMYDK6DCyzAXwPwHIDsRApvB3bOgjFORxdWdTsvlreUH+etwMgDB3i5MuL0Fd7auI91ThRAOScJIFaNYrWQsXgyA7+DBpX7/1sSt/Bb9Gzp0vNrpVQz6SljyXVHgt7EQtw1cfdUbVl3lp3N7yzdZmbH2OF+si8ZkseGk1/F49wie7dkYTxf5ry6EcFzyHaqKyd20CUtiInofH7xu6VWq95qtZiZtmQTAPU3vITIgsiJKvNSWL2D3D6DTw+C5svidnSmKwsoDSbz96yHiM/IB6NYogDf6RdIoyFPj6oQQ4tokjFQxGT8vAMCnXz/0LqVr0/3doe+IzozGz9WPZ697tiLKu9SJv2HluXtDek+Chj0r57w1xImUHN5YdoD1x1IBCPV14393NJcF7YQQVYqEkSrEkppK9t9/A+A7qHS9RZJyk/hizxcAjIkag49LJfSVSDsBCx4BxQpthkKnpyr+nDVEbqGFz9YcZ86GaMxWBWeDnidvaMDIGxvh5lwJl96EEMKOJIxUIZlLl4LFgmub1rg2bVKq907ZNoV8Sz7XBV1Hv4b9KqjCixRkwU9DoSADQtvDHVNBflIvN0VRWL43kXd+O0RSVgEANzUN5PU7I6kf4KFxdUIIUTYSRqoIRVEuWhSvdKMiG+M3svr0agw6AxM6TkCvq+BZLDYbLH4CUg6DVx249wcwulbsOWuAI0nZvL5sP/9GpwMQ5ufG63dE0qtF2ZYCEEIIRyFhpIrI374d06lT6N3d8bntthK/r9BayDtb3gFgaPOhNPVrWlElXrD2HTj6OxhcYMgP4FW74s9ZjWUVmPnkz2N8vekUVpuCi5Oep29qxBM9GuBqlEsyQoiqT8JIFXG+46r37beh9yj5cPzc/XOJyY4h0C2QkW1GVlR5FxxYAuveVx/3+xTqRlX8OaspRVFYvDOeyb8fJjWnEIA+kcG8ensLwvzcNa5OCCHsR8JIFWDNzCTrj5VA6S7RxGbH8uW+LwF48foX8XSu4GmeiXvhl3OBp/Mz0KYSVwGuZg4kZPL60gNsP30WgIgAD97oF8kNTWRBQSFE9SNhpArI/PVXlMJCXJo0wbV1yRaVUxSFd7e+S6G1kI61O3Jr/VsrtsicFPWGVXMeNLwZbnmrYs9XTWXmmflw9RG+//c0NgXcnQ0827Mxj3Wrj4uTXJIRQlRPEkYc3H9vXC1p74i/Y/9mXdw6nPROvNLplYrtOWExwc8PQWYs+DWEQXOgMjq7ViM2m8LP22OZsvII6bkmAO5oXYcJtzenjo+bxtUJIUTFkjDi4AoOHKTw8GF0zs749LuzRO/Jt+Tz3tb3AHi4xcM08GlQkSXCHy9DzCZw8Yb7fgK3WhV7vmpmT2wGry3dz564TAAaB3nyZr9IujQK0LgyIYSoHBJGHFzGArXjqtctt2Dw9S3Re2bvnU1CbgJ1POrwROsnKrA6YNsc2P4VoIOBX0Jg6fqf1GTpuSbeX3mYn7bFoijg6eLE6F6NebhLfYwGWURQCFFzSBhxYLa8PLJ+/RUo+aJ4JzNPMvfAXABevv5l3I0VOOvi1Ab4/SX18c2vQZM+FXeuasRqU5i3NYYPVh4hM98MwN3XhTKubzOCvKUfixCi5pEw4sCy/liJLTcXY3g47h2uv+b+iqLwzpZ3sNgsdA/tTs/wClwHJiNGvU/EZoGWg6DbmIo7VzWy43Q6//vlAAcTswBoVtuLtwe05Pr6fhpXJoQQ2pEw4sDOX6LxHTgQnf7aw/YrT6/k38R/cdY7M77D+Iq7adWUCz8Ohbw0qNMG+n0mrd6vISW7kHd/P8yinXEAeLs68UKfpgztEI6TXJIRQtRwEkYcVOHx4+Tv2gUGAz53Dbjm/rnmXN7fqjYbG9ZqGGHeYRVTmKKovUTO7AOPQLh3HjhLA64rsVhtfLv5NFNXHyW70ALAkPZhvHhrUwI8S7fqshBCVFcSRhxUxsJFAHjeeCPGoKBr7j9j9wyS85MJ8wrjsVaPVVxh6z6Ag7+A3ghDvgefuhV3riru3+g0Xl96gCNnsgFoFerDW/0juS5cZhsJIcTFJIw4IJvJROYvvwDgO2jgNfc/dvYY3x/6HoDxHcbjYqign7gP/wZ/T1Qf3/4hhHeqmPNUccnZBUz89RDL9iQA4Otu5KU+zRhyfRgGvVzOEkKI/5Iw4oBy/voLa0YGTsHBeHbvfs39p+6YilWxcnP4zXSve+39yyT5kLoSL0CHJyDq4Yo5TxVmsyn8tC2Wyb8fIrvAgk4H93cM5/lbmlLLw1nr8oQQwmFJGHFA5zuu+tx9Fzqnq/8V7U7ezfr49Rh0BsZGja2YgvLS4cf7wJQD9btDn3cq5jxV2LEz2YxfvK9oLZlWoT68c1crWtX10bgyIYRwfBJGHIwpLo7cTZsAdRbNtXy++3MA+jXsR7h3uP0Lslpg4aNw9iT4hsPgb8BgtP95qqgCs5Xpfx9nxj8nMFsV3J0NvNC7KQ93qS+XZIQQooQkjDiYjEXqjaseXbrgXPfqN4duT9rOv4n/4qR34sk2T1ZMQatehei1YPRQW717+FfMeaqgzSfSmLBkH9GpuQDc3CyItwa0JNRX1pIRQojSkDDiQBSLhczFSwDwHTzo6vsqCtN2TwPg7kZ3E+oZav+Cdn0PW2aoj++eCcGR9j9HFZSRZ+KdFYf4ebvaMyTIy4U3+kXSt2Xtil2QUAghqikJIw4kZ/16LGfOYPD1xfPmm6+675akLew4swOj3sjw1sPtX0zsVvj1XFfVG8dD85It0ledKYrCsj0JvLX8IGnnVta9v2M4L93aDB83uXQlhBBlJWHEgZzvLeIzYAB65yvPvlAUhWm71FGRwU0GU9ujtn0LyUqA+Q+A1QTN7oAeL9n3+FVQTFoeE37Zx/pjqYC6su7ku1vRXtq4CyFEuUkYcRDm5GRy1q4Frt1bZEP8Bvak7MHF4MKwVsPsXEg+/DQUcs5AUCTcNRNK0Iq+ujJbbXy14SRT/zxKgdmGs5Oe53o24okeDXF2qrl/LkIIYU8SRhxE5i9LwWrF7brrcGnU6Ir7KYpSNIPm3qb3EugeaL8iFAWWj4KEXeDmB/fNAxdP+x2/itkTm8G4xfs4dG5Ru84N/Jl0V0saBNbcPxMhhKgIEkYcgGKzkbFQ7S3iO+jqN67+Hfs3B9IO4ObkZv+275s+g73zQWeAe76BWvXte/wqIqfQwgcrj/DN5lMoitpBdcJtzRkUVVduUBVCiAogYcQB5G3dhjkmBr2HB959b73ifjbFxvTd0wG4v/n9+Lna8X6FY3/Cn6+rj/u+BxE97HfsKmT1wTO8tnQ/iZkFANx1XSiv3t4cf1nUTgghKoyEEQdwflTE+4470LtfeQXcP0//yZGzR/AwevBwCzu2Y089BgsfA8UG7R6C6+18H0oVcCargNeXHuCPA0kAhPu5M+mulnRvbMfLYEIIIS5LwojGrBkZZK9aBVz9Eo3VZi0aFXmwxYP4uvrap4CCTLXVe2EmhHWC2z6EGnQpwmZT+GFrDFN+P0x2oQWDXscTPRrwXM/GuDkbtC5PCCFqBAkjGstcthzFZMKleXNcW165qdgfp/7gROYJvJy9eLDFg/Y5uc0Ki4ZB2jHwrgtDvgOnmrOg25GkbMYv3svOmAwA2oT5MvmuVrQI8da2MCGEqGEkjGhIURQyFiwA1Om8V7o50mKzMGOP2gn1kchH8Ha204flX2/BsVXg5Ar3/gCeQfY5roMrMFv5bM0xZv4TjcWm4OFs4KVbm/FAp3qynowQQmhAwoiGCvbupfDYMXQuLvjccccV9/s1+ldOZ53G18WX+5vfb5+T710AGz9WH/f/HELa2ue4Dm7j8VQmLNnHqbQ8AHq3CObN/pHU8ZH1ZIQQQisSRjRUdOPqrX0w+Fx+qXmzzcwXe74A4LGWj+Fh9Cj/ieN3wrJn1MfdxkCrq08nrg7Sc01M+u0Qi3aq68kEe7vwZr+W3NrSzt1rhRBClJqEEY1Yc3LJ/G0FcPUbV5ceX0p8Tjx+rn4MaTqk/CfOPqO2ercUQOM+0PN/5T+mA1MUhSW74nn714OczTOj08GDnerxQp+meLvKejJCCOEIJIxoJOv3FSh5eTjXr49b+/aX3cdkNTFz70wAhrUahrvxytN+S8RSCD8/CFnxENAEBs4GffWdMXIqNZcJv+xj4/E0AJoGezF5YCvahdfSuDIhhBAXkzCikaKOq4MHXfHG1UXHFpGUm0SQWxD3NL2nfCdUFPjteYjdAq4+cN9P6u/VkNlqY9a6aD796xiFFhsuTnpG9WrM8O4NMBpkPRkhhHA0EkY0UHDkKAV79oKTEz79+19+H0sBs/fOBmB46+G4GMrZAXTrLNj1Hej0MOgr8G9YvuM5qJ0xZ3ll8T4OJ2UD0LWRP5MGtKJ+gB3utRFCCFEhJIxo4PyoiNdNN+EUEHDZfX4+8jMp+SnU8ajD3Y3vLt8Jo/+BP8arj295Cxr1Kt/xHFB2gZn3Vx7hu39PoyhQy93I/+5owV3Xhcp6MkII4eAkjFQyW2EhmcuWAeolmsvJM+cxZ/8cAJ5s/STOhnI0Iks/CQseBsUKre+Fzs+U/VgO6o/9Sby+bD9nsgoBGNiuLhNub46fR81p4CaEEFWZhJFKlr36T2yZmTjVqYNH166X3eenIz+RXpBOXc+69GvUr+wnK8yGn4ZC/lkIaQd3flKtWr0nZubz2tIDrD54BoD6/u5MuqsVXRtdfrRJCCGEY5IwUsmKOq7efTc6w6UzWXLNuczdPxeAEW1GYNSXcfqpzQZLRkDyQfCsrXZYNbqWuW5HYrUpfLf5FO+vPEKuyYqTXseTNzTg2Z6NcTVW39lBQghRXUkYqUSm06fJ27IFdDp8777rsvt8f/B7MgozqO9dn9sb3F72k/3zHhz+FQzOMOR78A4p+7EcyKHELMYt3see2AwA2oX7Mvnu1jSt7aVtYUIIIcpMwkglyli0GACPbt0whoZe8nqWKYtvDn4DwFNtnsJJX8a/noNL4Z931cd3fgJh15ftOA4k32Tlk7+OMXt9NFabgpeLEy/1bcb9HcLRy3oyQghRpZWp6cL06dOJiIjA1dWVqKgo1q9ff8V9Fy9ezC233EJgYCDe3t507tyZlStXlrngqkqxWMhYooaRK3Vc/fbAt2Sbsmnk24hbI24t24mS9quXZwA6PQ1th5btOA5k3dEUen/8D1/8cwKrTaFvy9r8+fwNPNipngQRIYSoBkodRubPn8/o0aOZMGECu3btonv37vTt25eYmJjL7r9u3TpuueUWVqxYwY4dO7jpppu488472bVrV7mLr0py/vkHa0oqBj8/vG668ZLXMwoy+P7Q9wCMbDsSva4MOTE3DX66D8x50OAmdRpvFZaaU8jon3bx0FdbiU3Pp46PK7Mfas+MB6II9q4e978IIYQAnaIoSmne0LFjR9q1a8eMGTOKtjVv3pwBAwYwefLkEh0jMjKSIUOG8Nprr5Vo/6ysLHx8fMjMzMTb27s05TqM2BFPkbN2LX6PP0bwiy9e8vrHOz5mzv45NK3VlJ/v/Ln0YcRqhu/uglProVYEDF8D7n52qr5yKYrCgh1xvLPiEBnn1pN5uHN9XujTFE8XubIohBBVRUk/v0v1nd1kMrFjxw7GjRtXbHvv3r3ZtGlTiY5hs9nIzs7Gz+/KH5SFhYUUFhYWPc/KyipNmQ7HfOYMOevWAeA78NJLNGn5acw7PA+Ap9s+XbZRkT/Gq0HE2Utt9V5Fg0h0Sg4Tluxnc7S6nkzzOt5MvrsVbcN8tS1MCCFEhSlVGElNTcVqtRIcHFxse3BwMElJSSU6xocffkhubi733HPltVYmT57Mm2++WZrSHFrm4sVgs+HWPgqXBhGXvP7V/q/It+TT0r8lN4bdWPoT7PwOts0GdOrid0HNyl1zZbNYbcxcF80nfx3DZLHhatQzulcTHu8WIevJCCFENVemMe//ttdWFKVELbd//PFH3njjDZYuXUpQUNAV9xs/fjxjx44tep6VlUVYWFhZStWcYrORsXARALUGD77k9eS8ZOYfmQ/A09c9XfrW5Un7YcUL6uOeE6Bp33LVq4WDCVm8tGgP++PVEbDujQOYNKAV4f7lXKVYCCFElVCqMBIQEIDBYLhkFCQ5OfmS0ZL/mj9/Po8//jgLFiygV6+rr43i4uKCi0s5F4ZzEHn//os5Ph69lxdevXtf8vqX+76k0FpI28C2dA25fEfWKyrMVlu9WwqgcW/o9rydqq4cJouNaX8fZ/rfx7HYFLxdnXj9zkjubifryQghRE1SqvFvZ2dnoqKiWL16dbHtq1evpkuXLld8348//sgjjzzCvHnzuP32cjTyqoLOnuu46nPnHejd3Iq9lpSbxMKj6qJ5z1z3TOk+gBUFlo+CtOPgXRfumgn6qnM5Y29cBnd+toFP/zqGxabQu0Uwf469gYFRdSWICCFEDVPqyzRjx47lwQcfpH379nTu3JlZs2YRExPDiBFqb4vx48cTHx/Pt99+C6hB5KGHHuKTTz6hU6dORaMqbm5u+Pj42PFLcTyWs2fJ/vMvAHwvc4lm1t5ZmG1m2ge3p0PtDqU7+PavYP8i0DvB4LlV5obVArOVj/88xqx1J7Ap4OfhzFv9I7m9VR0JIUIIUUOVOowMGTKEtLQ03nrrLRITE2nZsiUrVqygXr16ACQmJhbrOTJz5kwsFgtPP/00Tz/9dNH2hx9+mK+//rr8X4EDy1y6FMxmXCMjcW3evNhrcdlxLDm2BCjDqEjiHnX2DECvNyCslEFGI9tPpfPSor1Ep+QC0K9NCK/f2QJ/z+pxSU4IIUTZlOkG1pEjRzJy5MjLvvbfgLF27dqynKLKUxSFjAXqJRjfwZdO5525dyYWxULnOp2JCo4q+YELMuHnh8FaCE1vg87P2KvkCpNnsjDljyN8s/kUigJBXi5MHNCS3pG1tS5NCCGEA5AOUhUkf9duTCdOoHNzw/s/98mczjrN8hPLAXVUpMQUBZY+A2dPgk84DJgODn5pY9PxVF5evJfY9HwABkfV5dXbW+DjXsbViIUQQlQ7EkYqSMZCdVTE+9ZbMXgVX1F2xp4ZWBUrPer2oHVg65IfdOssOLQM9EYY/DW41bJjxfaVVWBm8orD/LhVvWQX4uPK5IGtuaFJoMaVCSGEcDQSRiqANSeHrN9/By69RBOdEc2K6BWA2m21xOJ3wMoJ6uPeE6FuKS7tVLK/jyTzyuJ9JGYWAPBAp3BevrUZXq4yGiKEEOJSEkYqQNavv6Hk5+PcsCFu111X7LXpe6ajoNAzrCct/FuU7ID5Z2HBI2AzQ/N+0PFJ+xdtBxl5Jt769SCLd8YDEO7nznsDW9O5ob/GlQkhhHBkEkYqwPlLNL6DBhWbJXMk/QgrT60E1JV5S0RR4JenISMGatWH/tMc8j6RP/Yn8eov+0nNKUSng8e6RvB87ya4O8s/MSGEEFcnnxR2VnDoEAX794PRiE//fsVem757OgB96vehqV/Tkh3w3+lw5DcwOMPgb8DVsXqzpOYU8vqyA/y2NxGAhoEeTBnUhqh6jns/ixBCCMciYcTOzk/n9br5ZpwuWpn4QNoB1sSuQa/TM7JNCUdFYrfB6tfUx33egZC2dq627BRFYfneRN5YdoD0XBMGvY4nezTguZsb42o0aF2eEEKIKkTCiB3ZCgrIXK5O2f3vjauf7/ocgNsibqOBb4NrHywvHRY+CjYLRN4N1w+ze71ldSargFd/2c/qg2cAaFbbi/cHtaFVXccatRFCCFE1SBixo+xVq7BlZ2MMCcGjc+ei7XtS9rA+fj0GnYERbUZc+0A2GywZAZmx4NcQ7vzEIe4TURSFhTviePvXg2QVWDAadDxzU2OeurEhzk5VZ10cIYQQjkXCiB1l/HxuUbxBA9FdtGjd+VGROxveST3vetc+0KZP4dhKMLio/URcvSui3FKJO5vHK0v2s+5oCgCt6/owZVBrmtXWvjYhhBBVm4QROyk8eZK87dtBr8f37ruLtu84s4PNiZtx0jnxZOsSTMk9vRn+ekt93Pc9qFOKpmgVwGZT+GFrDO+uOESuyYqzk56xtzRhWLcInAwyGiKEEKL8JIzYSeaiRQB4du+Osba65oqiKEzbNQ2AuxrfRV2vulc/SG4qLHwMFCu0ugeiHqnIkq/pVGouLy/ay5aT6QBE1avFlEGtaRjoqWldQgghqhcJI3agmExkLPkFKH7j6pakLWw/sx2j3sgTrZ+4+kFsNlj8BGQnQEATuGOqZveJWG0Kczee5INVRygw23AzGnjp1qY81Lk+Br32964IIYSoXiSM2EH22rVY09IwBAbgecMNQPFRkcFNBlPb4xor1G74CE78BU5uaj8RF21GH44nZ/Piwr3siskAoHMDf94b2Jpwf3dN6hFCCFH9SRixg6KOqwPuQmdU11/ZmLCRPSl7cDG4MKzVNablntoAf09SH9/+AQSXsE28HVmsNmaui+aTP49hstrwdHHilduac1+HsGJdZIUQQgh7kzBSTuaEBHLXbwDAd6B64+rFoyJDmg4h0P0qK9XmJMPCx0GxQZuhcN0DFV7zfx1MyOKlRXvYH58FwI1NA3nnrlaE+LpVei1CCCFqHgkj5ZSxeAkoCu4dOuBcvz4Aa2PXciDtAG5ObjzW8rErv9lmhcXDIScJApupoyKVyGSxMe3v40z/+zgWm4K3qxOv3xnJ3e1CZTRECCFEpZEwUg6K1UrGYnUWje/gwQDYFBuf71b7igxtNhR/t6usWLvufYheC0Z39T4RZ4+KLrnI3rgMXlywlyNnsgHo3SKYiQNaEuTtWmk1CCGEECBhpFxyN23CkpCI3scHr963APDn6T85cvYIHkYPHol85Mpvjl4La99VH98xFYKaVXi9AAVmKx//eYxZ605gU8DPw5m3+kdye6s6MhoihBBCExJGyuH8ong+d96J3sUFq81atDLvgy0exNfV9/JvzE6CRcMABa57ENrcWyn1bj+VzkuL9hKdkgtAvzYhvH5nC/w9XSrl/EIIIcTlSBgpI0taGtlr1gAXeousPLWSE5kn8HL24sEWD17+jVaLGkRyUyAoEm57v8JrzTNZmPLHEb7ZfApFgSAvFyYOaEnvyGtMNxZCCCEqgYSRMsr85RewWHBt3RrXpk2x2CzM2DMDgIdbPIy38xXWbPnnXTi1Hpw94Z5vwFixM1Y2HU/l5cV7iU3PB2BwVF1evb0FPu7GCj2vEEIIUVISRspAUZSiSzS+gwYC8Fv0b5zKOoWviy8PtLjC9Nzjf8G6czNm7vwEAhpXWI1ZBWYmrzjMj1tjAAjxcWXywNbc0OQq04yFEEIIDUgYKYP8HTswnTqFzt0d79tux2wz88WeLwB4tOWjeBgvMysmK0GdxosC7R+DVoMu3cdO/j6SzCuL95GYWQDAA53CefnWZni5ymiIEEIIxyNhpAzOj4p439YXg6cHS44uJC4nDj9XP+5tepmbUa0WtbFZXhrUbg19JldMXXkm3vr1IIt3xgMQ7ufOewNb07nhVaYXCyGEEBqTMFJK1qwsslauBKDWoEGYrCZm7p0JwLBWw3A3XmYNl78nQswmcPaCwV+D0f69PP7Yn8Srv+wnNacQnQ4e6xrB872b4O4sf8VCCCEcm3xSlVLmr7+iFBTg0rgxrm3aMP/IfJJykwhyC2Jwk8GXvuHoStgwVX3cfxr4N7RrPem5Jl5fdoDlexIAaBjowZRBbYiqV8uu5xFCCCEqioSRUipaFG/wIAqthczeOxuAYa2H4er0nxGPjFhY8qT6uMMTEDnArrWooyH7SM0xYdDreKJHA0bd3BhXo8Gu5xFCCCEqkoSRUsjff4DCg4fQGY1433knPx5dQHJ+MrU9ajOw8cDiO1vNsPAxyD8LIddB74l2q+Nsrok3lh9g6W51NKRxkCcfDG5DmzBfu51DCCGEqCwSRkohY+ECALxuuQWTpwtf7vsSgCdbP4mzwbn4zn++AXFbwcVHvU/EyT5dTlcdSOKVJeq9IXodPHlDQxkNEUIIUaVJGCkhW14eWb/+BoDvPYP56chPpBekU9ezLv0b9S++8+EVsHma+njA51CrfrnPn5Fn4s3lB1myS50p0zDQgw/vaUtbGQ0RQghRxUkYKaGsP1Ziy8nBGBaGcl0kc5e8CMCINiMw6i/q33H2NPwyQn3c6Wlofme5z/3nwTOMX7KPlGx1NGR4jwaM6dVERkOEEEJUCxJGSqjoxtWBA5l35EcyCjOo712f2xvcfmEniwkWPgoFmRDaHnq9Ua5zZuaZefPXA0V9QxoEevDB4Da0C5eZMkIIIaoPCSMlUHjiBPk7d4LBgOGOXny9Xl0Eb0SbETjpL/ojXP0axO8AV18YPBecnC9/wBJYc/gM4xfv40yW2jdkePcGjL1FRkOEEEJUPxJGSuB8x1XPG25gXuofZJuyaejTkFvr33php4PLYIu6UB53zQTf8DKdKzPfzNu/HmThjjgAGgR48P7g1kTV8yvX1yCEEEI4Kgkj12AzmchcuhQAY/9b+e7gJABGth2JQX9ulCI9GpY+rT7u8hw0vfVyh7qmv48kM37RPpKyCtDp4PGuEbzQp6mMhgghhKjWJIxcQ86aNVjPnsUpKIj5tY6Rm5RL01pN6VWvl7qDuQAWPAKFWRDWEW5+rdTnyCowM/HXg/y8XR0Nqe/vzvuD23B9fRkNEUIIUf1JGLmGjJ/V3iLO/W7lh2M/AfB026fR6/TqDqsmQOIecPODQXPBULqVcdcdTeHlRXtJzFRHQx7tEsGLfZri5iyjIUIIIWoGCSNXYYqLI3fTJgB+bZ5Pfmo+kf6R3Bh2o7rD/kWwTW18xt2zwCe0xMfOLjAz6bdD/LQtFoB6/u68P6gNHSJkNEQIIUTNImHkKjIXLwbA2LE9c8+uANRREZ1OB2knYNkodcduY6HxLSU+7vpjKby8cC8JmQUAPNKlPi/d2lRW2BVCCFEjyaffFSgWCxmL1DCyKcqNQmshbQLb0C20G5jz4eeHwZQN9brCTRNKdMycQguTfjvEj1tjAAj3c2fKoNZ0auBfYV+HEEII4egkjFxBzoYNWM6cQefjzTSvbQA8c90z6qjIH+PgzD5wD4CBc8Bw7T/GjcdTeWnhXuIz8gF4uHM9Xu7bTEZDhBBC1HjySXgF5zuuHulYh3z9CdoHt6dj7Y6w92fY8TWgg4Ffgnedqx4nt9DC5N8P8f2/6mhI3VpuTBnUmi4NAyr4KxBCCCGqBgkjl2FJSSHn77UAzK53Ejh3r0jqMVg+Wt3phpeg4U1XPc6mE+poSNxZdTTkwU71GNe3GR4u8scuhBBCnCefipeRseQXsFpJaejP6YBMOtXpRHu/FvDlzWDOhYgecMPLV3x/bqGF9/44zLebTwMQ6uvG+4Na06WRjIYIIYQQ/yVh5D8URSm6RLO4WSag3ivCihch+SB4BMHdX4L+8n1A/o1O48WFe4hNV0dD7u8YzvjbmuMpoyFCCCHEZckn5H/kbd2GOSYGk6sTG5opdA/tQZv4A7D7e9DpYdAc8Aq+9H0mC1P+OMLXm04B6mjIewNb062xjIYIIYQQVyNh5D8yFqgdV9c1s1HorOfperfD/MfVF28cr16i+Y+tJ9N5ceEeTqflAXBfhzBeua05Xq6l68YqhBBC1EQSRi5izcgge9UqAP5qo6NnaA8iV74JlnxocBN0f77Y/vkmK1NWHubrTadQFAjxceXdga3p0SRQi/KFEEKIKknCyEUyl/+KYjJxKghO1IH3zmZC6hHwqgN3zy52n8i2U+m8uGAPp86NhgxpH8aEO5rjLaMhQgghRKlIGDlHUZSiSzR/tdHT27cZTXcvV+8TGTgHPNXRjnyTlQ9WHeGrjSdRFKjt7cq7A1txY9MgLcsXQgghqiwJI+cU7NtH4dGjmJxgY6Se7w+rC+TR839QvysAO06n88KCvZxMzQVgcFRdXr2jBT5uMhoihBBClJWEkXMyFqjTef9tquMGg56GBbnQuDd0HU2B2cpHq48ye300igLB3i68e3drbmomoyFCCCFEeUkYAWy5uZz9dTk6YG1bPe+eiQPvUBjwBTvjMnlhwR6iU9TRkIHt6vLaHS3wcZfRECGEEMIeJIwAWb//ji6/gIRa0Mg3l3pnofCuL/loXTKz10VjUyDIy4V3B7aiZ7NLe4wIIYQQouwkjABx877BAKxto+PJzEzi27/Mw0vMHE+OBuDu60J5/c5IGQ0RQgghKoC+LG+aPn06ERERuLq6EhUVxfr166+6/z///ENUVBSurq40aNCAL774okzFVoSCo0cxHDyORQ+e9fPI9+pC9/WRHE/OIdDLhS8fas9HQ9pKEBFCCCEqSKnDyPz58xk9ejQTJkxg165ddO/enb59+xITE3PZ/U+ePMltt91G9+7d2bVrF6+88grPPfccixYtKnfx9nD462kA7GoEfS1GBiU+iE3RMaBtCKvH9KBXC7ksI4QQQlQknaIoSmne0LFjR9q1a8eMGTOKtjVv3pwBAwYwefLkS/Z/+eWXWbZsGYcOHSraNmLECPbs2cPmzZtLdM6srCx8fHzIzMzE29u7NOVelbWggN2do3DPt7HxdhO/Oj9HnEck79zVkt6Rte12HiGEEKImKunnd6lGRkwmEzt27KB3797Ftvfu3ZtNmzZd9j2bN2++ZP8+ffqwfft2zGbzZd9TWFhIVlZWsV8V4c9p43DPt5HmBdFuNxHeugerx/SQICKEEEJUolKFkdTUVKxWK8HBxS9dBAcHk5SUdNn3JCUlXXZ/i8VCamrqZd8zefJkfHx8in6FhYWVpswSsVktmFeq69BEN3Wmx5C3+PS+66jl4Wz3cwkhhBDiysp0A6tOpyv2XFGUS7Zda//LbT9v/PjxZGZmFv2KjY0tS5lXpTc44fHgoxxs5kTXsdO4tVWI3c8hhBBCiGsr1dTegIAADAbDJaMgycnJl4x+nFe7du3L7u/k5IS/v/9l3+Pi4oKLi0tpSiuTmx56ER56scLPI4QQQogrK9XIiLOzM1FRUaxevbrY9tWrV9OlS5fLvqdz586X7L9q1Srat2+P0SjTZYUQQoiartSXacaOHcuXX37JV199xaFDhxgzZgwxMTGMGDECUC+xPPTQQ0X7jxgxgtOnTzN27FgOHTrEV199xZw5c3jhhRfs91UIIYQQosoqdQfWIUP+397dhjTV/2EAv5ybj6TQk00npqGZRVaKphJCmEGBBEVCERYFSURmVBhFZgRRUZChBWL2RkuSgl5Y6ZtkWhRaQrTAUItELTTElT2p3/vF/V//e2l1n93bznZ2fWAvOv4m17hc+/rz7Cwfw8PDOHnyJAYGBrBkyRI0NjYiJiYGADAwMGB3zZHY2Fg0NjaiuLgYFRUViIyMRHl5OTZu3Oi8R0FEREReS/F1RtTgquuMEBERkeu45DojRERERM7GYYSIiIhUxWGEiIiIVMVhhIiIiFTFYYSIiIhUxWGEiIiIVMVhhIiIiFTFYYSIiIhUxWGEiIiIVKX4cvBqsF0kdnR0VOUkRERE9G/ZXrf/dLF3rxhGrFYrACA6OlrlJERERKSU1WpFeHj4L7/uFZ9NMzk5if7+fsyYMQN+fn5O+76jo6OIjo7G27dv+Zk3HoKdeBb24VnYh2dhH38mIrBarYiMjIRO9+szQ7xiZ0Sn08FkMrns+4eFhfEHycOwE8/CPjwL+/As7OP3frcjYsMTWImIiEhVHEaIiIhIVT49jAQGBqK0tBSBgYFqR6H/YSeehX14FvbhWdiH83jFCaxERESkXT69M0JERETq4zBCREREquIwQkRERKriMEJERESq0vwwUllZidjYWAQFBSElJQVms/m361taWpCSkoKgoCDExcXhypUrbkrqG5T0cevWLaxZswZz5sxBWFgYMjIycP/+fTem9Q1KnyM2bW1t0Ov1WLZsmWsD+hilfXz9+hVHjx5FTEwMAgMDsWDBAly9etVNabVPaR+1tbVITk5GSEgIjEYjduzYgeHhYTel9WKiYTdu3BCDwSBVVVVisVikqKhIQkND5c2bN9Ou7+npkZCQECkqKhKLxSJVVVViMBikoaHBzcm1SWkfRUVFcubMGXny5Il0dXXJkSNHxGAwyNOnT92cXLuUdmIzMjIicXFxkpubK8nJye4J6wMc6SMvL0/S09OlublZent75fHjx9LW1ubG1NqltA+z2Sw6nU4uXrwoPT09YjabZfHixbJhwwY3J/c+mh5G0tLSpLCw0O5YYmKilJSUTLv+8OHDkpiYaHds9+7dsnLlSpdl9CVK+5hOUlKSlJWVOTuaz3K0k/z8fDl27JiUlpZyGHEipX3cvXtXwsPDZXh42B3xfI7SPs6dOydxcXF2x8rLy8VkMrkso1Zo9s803759Q0dHB3Jzc+2O5+bm4uHDh9Pe59GjR1PWr127Fu3t7fj+/bvLsvoCR/r42eTkJKxWK2bOnOmKiD7H0U5qamrQ3d2N0tJSV0f0KY70cefOHaSmpuLs2bOIiopCQkICDh48iM+fP7sjsqY50kdmZib6+vrQ2NgIEcG7d+/Q0NCA9evXuyOyV/OKD8pzxNDQECYmJhAREWF3PCIiAoODg9PeZ3BwcNr14+PjGBoagtFodFlerXOkj5+dP38enz59wubNm10R0ec40smrV69QUlICs9kMvV6z/32owpE+enp60NraiqCgINy+fRtDQ0PYs2cPPnz4wPNG/iNH+sjMzERtbS3y8/Px5csXjI+PIy8vD5cuXXJHZK+m2Z0RGz8/P7t/i8iUY39aP91xcozSPmyuX7+OEydOoL6+HnPnznVVPJ/0bzuZmJjAli1bUFZWhoSEBHfF8zlKniOTk5Pw8/NDbW0t0tLSsG7dOly4cAHXrl3j7oiTKOnDYrFg3759OH78ODo6OnDv3j309vaisLDQHVG9mmZ/tZk9ezb8/f2nTLDv37+fMunazJs3b9r1er0es2bNcllWX+BIHzb19fXYuXMnbt68iZycHFfG9ClKO7FarWhvb8ezZ8+wd+9eAH+/GIoI9Ho9mpqasHr1ardk1yJHniNGoxFRUVF2H9G+aNEiiAj6+voQHx/v0sxa5kgfp0+fRlZWFg4dOgQAWLp0KUJDQ7Fq1SqcOnWKu+u/odmdkYCAAKSkpKC5udnueHNzMzIzM6e9T0ZGxpT1TU1NSE1NhcFgcFlWX+BIH8DfOyLbt29HXV0d/+7qZEo7CQsLw/Pnz9HZ2fnjVlhYiIULF6KzsxPp6enuiq5JjjxHsrKy0N/fj48fP/441tXVBZ1OB5PJ5NK8WudIH2NjY9Dp7F9W/f39Afx/l51+Qa0zZ93B9ras6upqsVgssn//fgkNDZXXr1+LiEhJSYls27btx3rbW3uLi4vFYrFIdXU139rrREr7qKurE71eLxUVFTIwMPDjNjIyotZD0BylnfyM76ZxLqV9WK1WMZlMsmnTJnnx4oW0tLRIfHy87Nq1S62HoClK+6ipqRG9Xi+VlZXS3d0tra2tkpqaKmlpaWo9BK+h6WFERKSiokJiYmIkICBAVqxYIS0tLT++VlBQINnZ2XbrHzx4IMuXL5eAgACZP3++XL582c2JtU1JH9nZ2QJgyq2goMD9wTVM6XPknziMOJ/SPl6+fCk5OTkSHBwsJpNJDhw4IGNjY25OrV1K+ygvL5ekpCQJDg4Wo9EoW7dulb6+Pjen9j5+Itw7IiIiIvVo9pwRIiIi8g4cRoiIiEhVHEaIiIhIVRxGiIiISFUcRoiIiEhVHEaIiIhIVRxGiIiISFUcRoiIiEhVHEaIiIhIVRxGiIiISFUcRoiIiEhVHEaIiIhIVX8BX/MArjyXKEoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "x=np.arange(10)*0.1\n",
    "for ind in [2,4,6,8]:\n",
    "    y=(1 - np.exp(-x* ind)) / (1 + np.exp(-x * ind))\n",
    "    plt.plot(x,y,label=str(ind))\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a6f4673",
   "metadata": {},
   "outputs": [],
   "source": [
    "class adjust_r():\n",
    "    def __init__(self,jc_scare,h,q):\n",
    "        self.jc_scare=jc_scare\n",
    "        self.h = h\n",
    "        self.q = q\n",
    "        self.env=environment()\n",
    "\n",
    "    def down_adjust(self, cd8_target, r_origin, error_origin,mu,j1,R_L_prop,a,ar_ratio,tumor,t1):\n",
    "        r_last = r_origin\n",
    "        error_last = error_origin\n",
    "        error = error_origin\n",
    "        while error <= error_last:\n",
    "            r = r_last - 0.002\n",
    "            jc_ = (1 - np.exp(-j1 * self.jc_scare)) / (1 + np.exp(-j1 * self.jc_scare))\n",
    "            x_ =self.env.fit_pre_treatment(mu, r, jc_, R_L_prop,  self.h,self.q,a,ar_ratio,t1)\n",
    "            j = cal_time(x_, tumor)\n",
    "            cd8_pred = min(x_[j][2] / (x_[j][0] + x_[j][1]), 1)\n",
    "\n",
    "            error_last = error\n",
    "            error = np.abs(cd8_pred - cd8_target)\n",
    "#             print(error,cd8_pred,cd8_target)\n",
    "            r_last = r\n",
    "            print(error,error_last,cd8_pred,cd8_target,r)\n",
    "        return r + 0.002\n",
    "\n",
    "\n",
    "    def up_adjust(self, cd8_target, r_origin, error_origin,mu,j1,R_L_prop,a,ar_ratio,tumor,t1):\n",
    "        r_last = r_origin\n",
    "        error_last = error_origin\n",
    "        error = error_origin\n",
    "        while error <= error_last:\n",
    "            r = r_last + 0.002\n",
    "            jc_ = (1 - np.exp(-j1 * self.jc_scare)) / (1 + np.exp(-j1 * self.jc_scare))\n",
    "            x_ = self.env.fit_pre_treatment(mu, r, jc_, R_L_prop,  self.h,self.q,a,ar_ratio,t1)\n",
    "            j = cal_time(x_, tumor)\n",
    "            cd8_pred = min(x_[j][2] / (x_[j][0] + x_[j][1]), 1)\n",
    "            error_last = error\n",
    "            error = np.abs(cd8_pred - cd8_target)\n",
    "            r_last = r\n",
    "            print(error,cd8_pred,cd8_target,r)\n",
    "            if r > 0.25:\n",
    "                break\n",
    "        return r - 0.002"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64d9b9df",
   "metadata": {},
   "outputs": [],
   "source": [
    "# given fixed j1_scare, q, h, estimate the latent parameter r\n",
    "r_l=np.ones(num_cases)*0.1\n",
    "tcga_all['time']=int(600)\n",
    "mae, cd8_l,time_l = opt_fix_para(r_l,6,0.01,0.03,num_cases,tcd8_target)\n",
    "\n",
    "r_estimator=adjust_r(6, 0.01, 0.03)\n",
    "for i in trange(len(cd8_l)):\n",
    "    tcd8_target_s = tcd8_target[i]\n",
    "    r_origin = r_l[i]\n",
    "    error_origin = np.abs(tcd8_target_s - cd8_l[i])\n",
    "    mu = tcga_all.loc[i, 'mu']\n",
    "    jc=tcga_all.loc[i,'mhc_norm']\n",
    "    R_L_prop=tcga_all.loc[i,'Treg_pro']\n",
    "    tumor=tcga_all.loc[i,'tumor_size']\n",
    "    a=tcga_all.loc[i,'a']\n",
    "    ar_ratio=tcga_all.loc[i,'ar_ratio']\n",
    "    t1=int(tcga_all.loc[i,'time'])\n",
    "    if tcd8_target_s > cd8_l[i]:\n",
    "        r = r_estimator.up_adjust(tcd8_target[i], r_origin, error_origin,mu,jc,R_L_prop,a,ar_ratio,tumor,t1)\n",
    "    else:\n",
    "        r = r_estimator.down_adjust(tcd8_target[i], r_origin, error_origin,mu,jc,R_L_prop,a,ar_ratio,tumor,t1)\n",
    "    r_l[i] = r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed4fac3e",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "mae, cd8_l,time_l= opt_fix_para(r_l,6, 0.01, 0.03,num_cases,tcd8_target)\n",
    "mae"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b782313",
   "metadata": {},
   "outputs": [],
   "source": [
    "# second round, the minimum mae 0.00514 with (j1_scale,h,q)= 6 0.01 0.03\n",
    "for j1_scale,h,q in test_range:\n",
    "    mae, cd8_l,time_l = opt_fix_para(r_l,j1_scale,h,q,num_cases,tcd8_target)\n",
    "    tcga_all['time']=time_l\n",
    "    print(j1_scale,h,q,mae)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "id": "6f5a99ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "#update r with j2=0.11\n",
    "env.update_fitted_j2(0.11)\n",
    "r_estimator.env.update_fitted_para(0.01,0.03,0.11)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "id": "c8815024",
   "metadata": {},
   "outputs": [],
   "source": [
    "r_l=tcga_all['r']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "id": "28829ab8",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████| 414/414 [03:02<00:00,  2.27it/s]\n"
     ]
    }
   ],
   "source": [
    "mae, cd8_l,time_l= opt_fix_para(r_l,6, 0.01, 0.03,num_cases,tcd8_target)\n",
    "mae"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3abfcb93",
   "metadata": {},
   "outputs": [],
   "source": [
    "r_estimator=adjust_r(6, 0.01, 0.03)\n",
    "for i in trange(len(cd8_l)):\n",
    "    tcd8_target_s = tcd8_target[i]\n",
    "    r_origin = r_l[i]\n",
    "    error_origin = np.abs(tcd8_target_s - cd8_l[i])\n",
    "    mu = tcga_all.loc[i, 'mu']\n",
    "    jc=tcga_all.loc[i,'mhc_norm']\n",
    "    R_L_prop=tcga_all.loc[i,'Treg_pro']\n",
    "    tumor=tcga_all.loc[i,'tumor_size']\n",
    "    a=tcga_all.loc[i,'a']\n",
    "    ar_ratio=tcga_all.loc[i,'ar_ratio']\n",
    "    t1=int(tcga_all.loc[i,'time'])\n",
    "    if tcd8_target_s > cd8_l[i]:\n",
    "        r = r_estimator.up_adjust(tcd8_target[i], r_origin, error_origin,mu,jc,R_L_prop,a,ar_ratio,tumor,t1)\n",
    "    else:\n",
    "        r = r_estimator.down_adjust(tcd8_target[i], r_origin, error_origin,mu,jc,R_L_prop,a,ar_ratio,tumor,t1)\n",
    "    r_l[i] = r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 182,
   "id": "9d6f4f9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "tcga_all['r']=r_l\n",
    "tcga_all.to_csv('tcga_all.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 183,
   "id": "c0a43f4c",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████| 414/414 [03:05<00:00,  2.23it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0029889716206676813\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "mae, cd8_l,time_l= opt_fix_para(r_l,6, 0.01, 0.03,num_cases,tcd8_target)\n",
    "print(mae)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
